matrix invertLT( const matrix &L )
{
   int N = L.size();
   matrix M( N, vector<double>( N, 0 ) );

   for ( int i = 0; i < N; i++ ) 
   {
      if ( L[i][i] == 0.0 )
      {
         cout << "*** Singular matrix ***\n";
         return M;
      }
      M[i][i] = 1.0 / L[i][i];
      for ( int j = 0; j < i; j++ )
      {
         for ( int k = j; k < i; k++ ) M[i][j] += L[i][k] * M[k][j];
         M[i][j] = -M[i][j] / L[i][i];
      }
   }
   return M;
}