00001
00006 inline long dsymatrix::dsysv(dgematrix& mat)
00007 {
00008 #ifdef CPPL_VERBOSE
00009 std::cerr << "# [MARK] dsymatrix::dsysv(dgematrix&)"
00010 << std::endl;
00011 #endif//CPPL_VERBOSE
00012
00013 #ifdef CPPL_DEBUG
00014 if(N!=mat.N){
00015 std::cerr << "[ERROR] dsymatrix::dsysv(dgematrix&) " << std::endl
00016 << "These two matrices cannot be solved." << std::endl
00017 << "Your input was (" << N << "x" << N << ") and ("
00018 << mat.N << "x" << mat.N << ")." << std::endl;
00019 exit(1);
00020 }
00021 #endif//CPPL_DEBUG
00022
00023 char UPLO('L');
00024 long NRHS(mat.N), LDA(N), *IPIV(new long[N]), LDB(mat.N), LWORK(-1), INFO(1);
00025 double *WORK( new double[1] );
00026 dsysv_(UPLO, N, NRHS, Array, LDA, IPIV, mat.Array, LDB, WORK, LWORK, INFO);
00027
00028 INFO=1;
00029 LWORK = long(WORK[0]);
00030 delete [] WORK; WORK = new double[LWORK];
00031 dsysv_(UPLO, N, NRHS, Array, LDA, IPIV, mat.Array, LDB, WORK, LWORK, INFO);
00032 delete [] WORK; delete [] IPIV;
00033
00034 if(INFO!=0){
00035 std::cerr << "[WARNING] dsymatrix::dsysv(dsymatrix&) "
00036 << "Serious trouble happend. INFO = " << INFO << "."
00037 << std::endl;
00038 }
00039 return INFO;
00040 }
00041
00042
00047 inline long dsymatrix::dsysv(dcovector& vec)
00048 {
00049 #ifdef CPPL_VERBOSE
00050 std::cerr << "# [MARK] dsymatrix::dsysv(dcovector&)"
00051 << std::endl;
00052 #endif//CPPL_VERBOSE
00053
00054 #ifdef CPPL_DEBUG
00055 if(N!=vec.L){
00056 std::cerr << "[ERROR] dsymatrix::dsysv(dcovector&) " << std::endl
00057 << "These matrix and vector cannot be solved." << std::endl
00058 << "Your input was (" << N << "x" << N << ") and ("
00059 << vec.L << ")." << std::endl;
00060 exit(1);
00061 }
00062 #endif//CPPL_DEBUG
00063
00064 char UPLO('L');
00065 long NRHS(1), LDA(N), *IPIV(new long[N]), LDB(vec.L), LWORK(-1), INFO(1);
00066 double *WORK( new double[1] );
00067 dsysv_(UPLO, N, NRHS, Array, LDA, IPIV, vec.Array, LDB, WORK, LWORK, INFO);
00068
00069 INFO=1;
00070 LWORK = long(WORK[0]);
00071 delete WORK; WORK = new double[LWORK];
00072 dsysv_(UPLO, N, NRHS, Array, LDA, IPIV, vec.Array, LDB, WORK, LWORK, INFO);
00073 delete [] WORK; delete [] IPIV;
00074
00075 if(INFO!=0){
00076 std::cerr << "[WARNING] dsymatrix::dsysv(dcovector&) "
00077 << "Serious trouble happend. INFO = " << INFO << "."
00078 << std::endl;
00079 }
00080 return INFO;
00081 }
00082
00086
00087
00094 inline long dsymatrix::dsyev(std::vector<double>& w, const bool& jobz=0)
00095 {
00096 #ifdef CPPL_VERBOSE
00097 std::cerr << "# [MARK] dsymatrix::dsyev(std::vector<double>&, const bool&)"
00098 << std::endl;
00099 #endif//CPPL_VERBOSE
00100
00101 w.resize(N);
00102 char JOBZ, UPLO('L');
00103 if(jobz==0){ JOBZ='N'; } else{ JOBZ='V'; }
00104 long LDA(N), INFO(1), LWORK(-1);
00105 double *WORK(new double[1]);
00106 dsyev_(JOBZ, UPLO, N, Array, LDA, &w[0], WORK, LWORK, INFO);
00107
00108 INFO=1;
00109 LWORK = long(WORK[0]);
00110 delete [] WORK; WORK = new double[LWORK];
00111 dsyev_(JOBZ, UPLO, N, Array, LDA, &w[0], WORK, LWORK, INFO);
00112 delete [] WORK;
00113
00114 if(INFO!=0){
00115 std::cerr << "[WARNING] dsymatrix::dsyev(vector<double>&, const bool&)"
00116 << std::endl
00117 << "Serious trouble happend. INFO = " << INFO << "."
00118 << std::endl;
00119 }
00120 return INFO;
00121 }
00122
00123
00130 inline long dsymatrix::dsyev(std::vector<double>& w, std::vector<dcovector>& v)
00131 {
00132 #ifdef CPPL_VERBOSE
00133 std::cerr << "# [MARK] dsymatrix::dsyev(std::vector<double>&, std::vector<dcovector>&)"
00134 << std::endl;
00135 #endif//CPPL_VERBOSE
00136
00137 w.resize(N); v.resize(N);
00138 for(long i=0; i<N; i++){ v[i].resize(N); }
00139 char JOBZ('V'), UPLO('L');
00140 long LDA(N), INFO(1), LWORK(-1);
00141 double *WORK(new double[1]);
00142 dsyev_(JOBZ, UPLO, N, Array, LDA, &w[0], WORK, LWORK, INFO);
00143
00144 INFO=1;
00145 LWORK = long(WORK[0]);
00146 delete [] WORK; WORK = new double[LWORK];
00147 dsyev_(JOBZ, UPLO, N, Array, LDA, &w[0], WORK, LWORK, INFO);
00148 delete [] WORK;
00149
00151 for(long i=0; i<N; i++){ for(long j=0; j<N; j++){
00152 v[j](i) = Array[i+N*j];
00153 }}
00154
00155 if(INFO!=0){
00156 std::cerr << "[WARNING] dsymatrix::dsyev"
00157 << "(vector<double>&, vector<dcovector>&)" << std::endl
00158 << "Serious trouble happend. INFO = " << INFO << "."
00159 << std::endl;
00160 }
00161 return INFO;
00162 }
00163
00164
00171 inline long dsymatrix::dsyev(std::vector<double>& w, std::vector<drovector>& v)
00172 {
00173 #ifdef CPPL_VERBOSE
00174 std::cerr << "# [MARK] dsymatrix::dsyev(std::vector<double>&, std::vector<drovector>&)"
00175 << std::endl;
00176 #endif//CPPL_VERBOSE
00177
00178 w.resize(N); v.resize(N);
00179 for(long i=0; i<N; i++){ v[i].resize(N); }
00180 char JOBZ('V'), UPLO('L');
00181 long LDA(N), INFO(1), LWORK(-1);
00182 double *WORK(new double[1]);
00183 dsyev_(JOBZ, UPLO, N, Array, LDA, &w[0], WORK, LWORK, INFO);
00184
00185 INFO=1;
00186 LWORK = long(WORK[0]);
00187 delete [] WORK; WORK = new double[LWORK];
00188 dsyev_(JOBZ, UPLO, N, Array, LDA, &w[0], WORK, LWORK, INFO);
00189 delete [] WORK;
00190
00192 for(long i=0; i<N; i++){ for(long j=0; j<N; j++){
00193 v[j](i) = Array[i+N*j];
00194 }}
00195
00196 if(INFO!=0){
00197 std::cerr << "[WARNING] dsymatrix::dsyev"
00198 << "(vector<double>&, vector<dcovector>&)" << std::endl
00199 << "Serious trouble happend. INFO = " << INFO << "."
00200 << std::endl;
00201 }
00202 return INFO;
00203 }
00204
00208
00209
00214 inline long dsymatrix::dsygv(dsymatrix& matB, std::vector<double>& w)
00215 {
00216 #ifdef CPPL_VERBOSE
00217 std::cerr << "# [MARK] dsymatrix::dsygv(dsymatrix&, std::vector<double>&)"
00218 << std::endl;
00219 #endif//CPPL_VERBOSE
00220
00221 #ifdef CPPL_DEBUG
00222 if(matB.N!=N){
00223 std::cerr << "[ERROR] dsymatrix::dsygv"
00224 << "(dsymatrix&, vector<double>&) " << std::endl
00225 << "The matrix B is not a matrix "
00226 << "having the same size as \"this\" matrix." << std::endl
00227 << "The B matrix is (" << matB.N << "x" << matB.N << ")."
00228 << std::endl;
00229 exit(1);
00230 }
00231 #endif//CPPL_DEBUG
00232
00233 w.resize(N);
00234 char JOBZ('N'), UPLO('L');
00235 long ITYPE(1), LDA(N), LDB(N), LWORK(-1), INFO(1);
00236 double *WORK(new double[1]);
00237 dsygv_(ITYPE, JOBZ, UPLO, N, Array, LDA, matB.Array, LDB, &w[0],
00238 WORK, LWORK, INFO);
00239 INFO=1;
00240 LWORK = long(WORK[0]);
00241 delete [] WORK; WORK = new double[LWORK];
00242 dsygv_(ITYPE, JOBZ, UPLO, N, Array, LDA, matB.Array, LDB, &w[0],
00243 WORK, LWORK, INFO);
00244 delete [] WORK;
00245
00246 if(INFO!=0){
00247 std::cerr << "[WARNING] dsymatrix::dsygv"
00248 << "(dsymatrix&, vector<double>&)"
00249 << std::endl
00250 << "Serious trouble happend. INFO = " << INFO << "."
00251 << std::endl
00252 << "Make sure that B matrix is a symmetric-definite."
00253 << std::endl;
00254 }
00255 return INFO;
00256 }
00257
00258
00263 inline long dsymatrix::dsygv(dsymatrix& matB, std::vector<double>& w,
00264 std::vector<dcovector>& v)
00265 {
00266 #ifdef CPPL_VERBOSE
00267 std::cerr << "# [MARK] dsymatrix::dsygv(dsymatrix&, std::vector<double>&, std::vector<dcovector>&)"
00268 << std::endl;
00269 #endif//CPPL_VERBOSE
00270
00271 #ifdef CPPL_DEBUG
00272 if(matB.N!=N){
00273 std::cerr << "[ERROR] dsymatrix::dsygv"
00274 << "(dsymatrix&, vector<double>&) " << std::endl
00275 << "The matrix B is not a matrix "
00276 << "having the same size as \"this\" matrix." << std::endl
00277 << "The B matrix is (" << matB.N << "x" << matB.N << ")."
00278 << std::endl;
00279 exit(1);
00280 }
00281 #endif//CPPL_DEBUG
00282
00283 w.resize(N);
00284 v.resize(N);
00285 char JOBZ('V'), UPLO('L');
00286 long ITYPE(1), LDA(N), LDB(N), LWORK(-1), INFO(1);
00287 double *WORK(new double[1]);
00288 dsygv_(ITYPE, JOBZ, UPLO, N, Array, LDA, matB.Array, LDB, &w[0],
00289 WORK, LWORK, INFO);
00290 INFO=1;
00291 LWORK = long(WORK[0]);
00292 std::cout << " LWORK = " << LWORK <<std::endl;
00293 delete [] WORK; WORK = new double[LWORK];
00294 dsygv_(ITYPE, JOBZ, UPLO, N, Array, LDA, matB.Array, LDB, &w[0],
00295 WORK, LWORK, INFO);
00296 delete [] WORK;
00297
00299 for(int i=0; i<N; i++){
00300 v[i].resize(N);
00301 for(int j=0; j<N; j++){
00302 v[i](j) =Darray[i][j];
00303 }
00304 }
00305
00306 if(INFO!=0){
00307 std::cerr << "[WARNING] dsymatrix::dsygv"
00308 << "(dsymatrix&, vector<double>&)"
00309 << std::endl
00310 << "Serious trouble happend. INFO = " << INFO << "."
00311 << std::endl
00312 << "Make sure that B matrix is a symmetric-definite."
00313 << std::endl;
00314 }
00315 return INFO;
00316 }