CPPScaLapack: dgematrix-lapack.hpp Source File

00001 //=============================================================================
00002 /*! solve A*X=Y using dgesv\n
00003   The argument is dgematrix Y. Y is overwritten and become the solution X.
00004   A is also overwritten and become P*L*U. */
00005 inline long dgematrix::dgesv(dgematrix& mat)
00006 {
00007 #ifdef  CPPSL_DEBUG
00008   if(M!=N || M!=mat.M){
00009     std::cerr << "[ERROR] dgematrix::dgesv(dgematrix&) " << std::endl
00010               << "These two matrices cannot be solved." << std::endl
00011               << "Your input was (" << M << "x" << N << ") and ("
00012               << mat.M << "x" << mat.N << ")." << std::endl;
00013     exit(1);
00014   }
00015 #endif//CPPSL_DEBUG 
00016   long NRHS(mat.N), IA(1), JA(1), *IPIV(new long[Ml+mb]), IB(1), JB(1), INFO(1);
00017   pdgesv_( N, NRHS, Array, IA, JA, Desc, IPIV,
00018            mat.Array, IB, JB, mat.Desc, INFO );
00019   delete [] IPIV;
00020   
00021   if(INFO!=0){
00022     std::cerr << "[WARNING] dgematrix::dgesv(dgematrix&) "
00023               << "Serious trouble happend. INFO = " << INFO << "."
00024               << std::endl;
00025   }
00026   return INFO;
00027 }
00028 
00029 //=============================================================================
00030 /*! solve A*x=y using dgesv\n
00031   The argument is dcovector y. y is overwritten and become the solution x.
00032   A is also overwritten and become P*L*U. */
00033 inline long dgematrix::dgesv(dcovector& vec)
00034 {
00035 #ifdef  CPPSL_DEBUG
00036   if(m!=n || m!=vec.l){
00037     std::cerr << "[ERROR] dgematrix::dgesv(dcovector&) " << std::endl
00038               << "These matrix and vector cannot be solved." << std::endl
00039               << "Your input was (" << m << "x" << n << ") and ("
00040               << vec.l << ")." << std::endl;
00041     exit(1);
00042   }
00043 #endif//CPPSL_DEBUG 
00044   long NRHS(1), IA(1), JA(1), *IPIV(new long[Ml+mb]), IB(1), JB(1), INFO(1);
00045   pdgesv_(N, NRHS, Array, IA, JA, Desc, IPIV,
00046           vec.Array, IB, JB, vec.Desc, INFO);
00047   delete [] IPIV;
00048   
00049   if(INFO!=0){
00050     std::cerr << "[WARNING] dgematrix::dgesv(dcovector&) "
00051               << "Serious trouble happend. INFO = " << INFO << "."
00052               << std::endl;
00053   }
00054   return INFO;
00055 }