9 std::cerr <<
"# [MARK] dgematrix::dgesv(dgematrix&)"
15 std::cerr <<
"[ERROR] dgematrix::dgesv(dgematrix&) " << std::endl
16 <<
"These two matrices cannot be solved." << std::endl
17 <<
"Your input was (" << M <<
"x" << N <<
") and ("
18 << mat.M <<
"x" << mat.N <<
")." << std::endl;
23 long NRHS(mat.N), LDA(N), *IPIV(
new long[N]), LDB(mat.M), INFO(1);
24 dgesv_(N, NRHS,
Array, LDA, IPIV, mat.Array, LDB, INFO);
28 std::cerr <<
"[WARNING] dgematrix::dgesv(dgematrix&) "
29 <<
"Serious trouble happend. INFO = " << INFO <<
"."
43 std::cerr <<
"# [MARK] dgematrix::dgesv(dcovector&)"
49 std::cerr <<
"[ERROR] dgematrix::dgesv(dcovector&) " << std::endl
50 <<
"These matrix and vector cannot be solved." << std::endl
51 <<
"Your input was (" << M <<
"x" << N <<
") and ("
52 << vec.L <<
")." << std::endl;
57 long NRHS(1), LDA(N), *IPIV(
new long[N]), LDB(vec.L), INFO(1);
58 dgesv_(N, NRHS,
Array, LDA, IPIV, vec.Array, LDB, INFO);
62 std::cerr <<
"[WARNING] dgematrix::dgesv(dcovector&) "
63 <<
"Serious trouble happend. INFO = " << INFO <<
"."
79 std::cerr <<
"# [MARK] dgematrix::dgels(dgematrix&)"
85 std::cerr <<
"[ERROR] dgematrix::dgels(dgematrix&) " << std::endl
86 <<
"These two matrices cannot be solved." << std::endl
87 <<
"Your input was (" << M <<
"x" << N <<
") and ("
88 << mat.M <<
"x" << mat.N <<
")." << std::endl;
95 for(
long i=0;
i<mat.M;
i++){
for(
long j=0; j<mat.N; j++){
103 long NRHS(mat.N), LDA(M), LDB(mat.M),
104 LWORK(min(M,N)+max(min(M,N),NRHS)), INFO(1);
105 double *WORK(
new double[LWORK]);
106 dgels_(TRANS, M, N, NRHS,
Array, LDA, mat.Array, LDB, WORK, LWORK, INFO);
111 for(
long i=0;
i<tmp.M;
i++){
for(
long j=0; j<tmp.N; j++){
119 std::cerr <<
"[WARNING] dgematrix::dgels(dgematrix&) "
120 <<
"Serious trouble happend. INFO = " << INFO <<
"."
132 std::cerr <<
"# [MARK] dgematrix::dgels(dcovector&)"
138 std::cerr <<
"[ERROR] dgematrix::dgels(dcovector&) " << std::endl
139 <<
"These matrix and vector cannot be solved." << std::endl
140 <<
"Your input was (" << M <<
"x" << N <<
") and ("
141 << vec.L <<
")." << std::endl;
148 for(
long i=0;
i<vec.L;
i++){ tmp(
i)=vec(
i); }
154 long NRHS(1), LDA(M), LDB(vec.L),
155 LWORK(min(M,N)+max(min(M,N),NRHS)), INFO(1);
156 double *WORK(
new double[LWORK]);
157 dgels_(TRANS, M, N, NRHS,
Array, LDA, vec.Array, LDB, WORK, LWORK, INFO);
162 for(
long i=0;
i<tmp.L;
i++){ tmp(
i)=vec(
i); }
168 std::cerr <<
"[WARNING] dgematrix::dgels(dcovector&) "
169 <<
"Serious trouble happend. INFO = " << INFO <<
"."
185 std::cerr <<
"# [MARK] dgematrix::dgels(dgematrix&, drovector&)"
191 std::cerr <<
"[ERROR] dgematrix::dgels(dgematrix&, drovector&) "
193 <<
"These two matrices cannot be solved." << std::endl
194 <<
"Your input was (" << M <<
"x" << N <<
") and ("
195 << mat.M <<
"x" << mat.N <<
")." << std::endl;
204 for(
long i=0;
i<mat.M;
i++){
for(
long j=0; j<mat.N; j++){
212 long NRHS(mat.N), LDA(M), LDB(mat.M),
213 LWORK(min(M,N)+max(min(M,N),NRHS)), INFO(1);
214 double *WORK(
new double[LWORK]);
215 dgels_(TRANS, M, N, NRHS,
Array, LDA, mat.Array, LDB, WORK, LWORK, INFO);
219 for(
long i=0;
i<residual.L;
i++){
for(
long j=0; j<M-N; j++){
220 residual(
i) += std::pow(mat(N+j,
i), 2.0);
224 for(
long i=0;
i<tmp.M;
i++){
for(
long j=0; j<tmp.N; j++){
232 std::cerr <<
"[WARNING] dgematrix::dgels(dgematrix&, drovector&) "
233 <<
"Serious trouble happend. INFO = " << INFO <<
"."
249 std::cerr <<
"# [MARK] dgematrix::dgels(dcovector&, double&)"
255 std::cerr <<
"[ERROR] dgematrix::dgels(dcovector&, double&) " << std::endl
256 <<
"These matrix and vector cannot be solved." << std::endl
257 <<
"Your input was (" << M <<
"x" << N <<
") and ("
258 << vec.L <<
")." << std::endl;
267 for(
long i=0;
i<vec.L;
i++){ tmp(
i)=vec(
i); }
273 long NRHS(1), LDA(M), LDB(vec.L),
274 LWORK(min(M,N)+max(min(M,N),NRHS)), INFO(1);
275 double *WORK(
new double[LWORK]);
276 dgels_(TRANS, M, N, NRHS,
Array, LDA, vec.Array, LDB, WORK, LWORK, INFO);
280 for(
long i=0;
i<M-N;
i++){ residual+=std::pow(vec(N+
i),2.0); }
283 for(
long i=0;
i<tmp.L;
i++){ tmp(
i)=vec(
i); }
289 std::cerr <<
"[WARNING] dgematrix::dgels(dcovector&, double&) "
290 <<
"Serious trouble happend. INFO = " << INFO <<
"."
305 const double RCOND =-1. )
308 std::cerr <<
"# [MARK] dgematrix::dgelss(dcovector&, dcovector&, long&, const double)"
314 std::cerr <<
"[ERROR] dgematrix::dgelss"
315 <<
"(dcovector&, dcovector&, long, double) " << std::endl
316 <<
"These matrix and vector cannot be solved." << std::endl
317 <<
"Your input was (" << M <<
"x" << N <<
") and ("
318 << B.L <<
")." << std::endl;
325 for(
long i=0;
i<B.L;
i++){ tmp(
i)=B(
i); }
332 long NRHS(1), LDA(M), LDB(B.L),
333 LWORK(3*min(M,N)+max(max(2*min(M,N),max(M,N)), NRHS)), INFO(1);
334 double *WORK(
new double[LWORK]);
335 dgelss_(M, N, NRHS,
Array, LDA, B.Array, LDB, S.Array,
336 RCOND, RANK, WORK, LWORK, INFO);
340 std::cerr <<
"[WARNING] dgematrix::dgelss"
341 <<
"(dcovector&, docvector&, long, const double) "
342 <<
"Serious trouble happend. INFO = " << INFO <<
"."
353 const double RCOND =-1. )
356 std::cerr <<
"# [MARK] dgematrix::dgelss(dgematrix&, dcovector&, long& const double)"
362 std::cerr <<
"[ERROR] dgematrix::dgelss"
363 <<
"(dgematrix&, dcovector&, long&, const double) " << std::endl
364 <<
"These matrix and vector cannot be solved." << std::endl
365 <<
"Your input was (" << M <<
"x" << N <<
") and ("
366 << B.M <<
"x" << B.N <<
")." << std::endl;
373 for(
long i=0;
i<B.M;
i++){
374 for(
long j=0; j<B.N; j++){
384 long NRHS(B.N), LDA(M), LDB(B.M),
385 LWORK(3*min(M,N)+max(max(2*min(M,N),max(M,N)), NRHS)), INFO(1);
386 double *WORK(
new double[LWORK]);
387 dgelss_(M, N, NRHS,
Array, LDA, B.Array, LDB, S.Array,
388 RCOND, RANK, WORK, LWORK, INFO);
392 std::cerr <<
"[WARNING] dgematrix::dgelss"
393 <<
"(dgematrix&, docvector&, long, const double) "
394 <<
"Serious trouble happend. INFO = " << INFO <<
"."
410 inline long dgematrix::dgeev(std::vector<double>& wr, std::vector<double>& wi)
413 std::cerr <<
"# [MARK] dgematrix::dgeev(std::vector<double>&, std::vector<double>&)"
419 std::cerr <<
"[ERROR] dgematrix::dgeev"
420 <<
"(vector<double>&, vector<double>&) " << std::endl
421 <<
"This matrix is not a square matrix." << std::endl
422 <<
"This matrix is (" << M <<
"x" << N <<
")." << std::endl;
427 wr.resize(N); wi.resize(N);
428 char JOBVL(
'N'), JOBVR(
'N');
429 long LDA(N), LDVL(1), LDVR(1), LWORK(3*N), INFO(1);
430 double *VL(NULL), *VR(NULL), *WORK(
new double[LWORK]);
431 dgeev_(JOBVL, JOBVR, N,
Array, LDA, &wr[0], &wi[0],
432 VL, LDVL, VR, LDVR, WORK, LWORK, INFO);
433 delete [] WORK;
delete [] VL;
delete [] VL;
436 std::cerr <<
"[WARNING] dgematrix::dgeev"
437 <<
"(vector<double>&, vector<double>&)"
439 <<
"Serious trouble happend. INFO = " << INFO <<
"."
453 inline long dgematrix::dgeev(std::vector<double>& wr, std::vector<double>& wi,
454 std::vector<dcovector>& vrr,
455 std::vector<dcovector>& vri)
458 std::cerr <<
"# [MARK] dgematrix::dgeev(std::vector<double>&, std::vector<double>&, std::vector<dcovector>&, std::vector<dcovector>&)"
464 std::cerr <<
"[ERROR] dgematrix::dgeev"
465 <<
"(vector<double>&, vector<double>&, "
466 <<
"vector<dcovector>&, vector<dcovector>&) "
468 <<
"This matrix is not a square matrix." << std::endl
469 <<
"This matrix is (" << M <<
"x" << N <<
")." << std::endl;
474 wr.resize(N); wi.resize(N); vrr.resize(N); vri.resize(N);
475 for(
long i=0;
i<N;
i++){ vrr[
i].resize(N); vri[
i].resize(N); }
477 char JOBVL(
'N'), JOBVR(
'V');
478 long LDA(N), LDVL(1), LDVR(N), LWORK(4*N), INFO(1);
479 double *VL(NULL), *WORK(
new double[LWORK]);
480 dgeev_(JOBVL, JOBVR, N,
Array, LDA, &wr[0], &wi[0],
481 VL, LDVL, VR.Array, LDVR, WORK, LWORK, INFO);
482 delete [] WORK;
delete [] VL;
485 for(
long j=0; j<N; j++){
486 if(fabs(wi[j])<1e-10){
487 for(
long i=0;
i<N;
i++){
488 vrr[j](
i) = VR(
i,j); vri[j](
i) = 0.0;
492 for(
long i=0;
i<N;
i++){
493 vrr[j](
i) = VR(
i,j); vri[j](
i) = VR(
i,j+1);
494 vrr[j+1](
i) = VR(
i,j); vri[j+1](
i) =-VR(
i,j+1);
501 std::cerr <<
"[WARNING] dgematrix::dgeev"
502 <<
"(vector<double>&, vector<double>&, "
503 <<
"vector<dcovector>&, vector<dcovector>&) "
505 <<
"Serious trouble happend. INFO = " << INFO <<
"."
519 inline long dgematrix::dgeev(std::vector<double>& wr, std::vector<double>& wi,
520 std::vector<drovector>& vlr,
521 std::vector<drovector>& vli)
524 std::cerr <<
"# [MARK] dgematrix::dgeev(std::vector<double>&, std::vector<double>&, std::vector<drovector>&, std::vector<drovector>&)"
530 std::cerr <<
"[ERROR] dgematrix::dgeev"
531 <<
"(vector<double>&, vector<double>&, "
532 <<
"vector<drovector>&, vector<drovector>&) "
534 <<
"This matrix is not a square matrix." << std::endl
535 <<
"This matrix is (" << M <<
"x" << N <<
")." << std::endl;
540 wr.resize(N); wi.resize(N); vlr.resize(N); vli.resize(N);
541 for(
long i=0;
i<N;
i++){ vlr[
i].resize(N); vli[
i].resize(N); }
543 char JOBVL(
'V'), JOBVR(
'N');
544 long LDA(N), LDVL(N), LDVR(1), LWORK(4*N), INFO(1);
545 double *VR(NULL), *WORK(
new double[LWORK]);
546 dgeev_(JOBVL, JOBVR, N,
Array, LDA, &wr[0], &wi[0],
547 VL.Array, LDVL, VR, LDVR, WORK, LWORK, INFO);
548 delete [] WORK;
delete [] VR;
551 for(
long j=0; j<N; j++){
552 if(fabs(wi[j])<1e-10){
553 for(
long i=0;
i<N;
i++){
554 vlr[j](
i) = VL(
i,j); vli[j](
i) = 0.0;
558 for(
long i=0;
i<N;
i++){
559 vlr[j](
i) = VL(
i,j); vli[j](
i) =-VL(
i,j+1);
560 vlr[j+1](
i) = VL(
i,j); vli[j+1](
i) = VL(
i,j+1);
568 std::cerr <<
"[WARNING] dgematrix::dgeev"
569 <<
"(vector<double>&, vector<double>&, "
570 <<
"vector<drovector>&, vector<drovector>&) "
572 <<
"Serious trouble happend. INFO = " << INFO <<
"."
592 std::vector<double>& wr, std::vector<double>& wi)
595 std::cerr <<
"# [MARK] dgematrix::dggev(dgematrix&, std::vector<double>&, std::vector<double>&)"
601 std::cerr <<
"[ERROR] dgematrix::dggev"
602 <<
"(dgematrix&, vector<double>&, vector<double>&) " << std::endl
603 <<
"This matrix is not a square matrix." << std::endl
604 <<
"This matrix is (" << M <<
"x" << N <<
")." << std::endl;
607 if(matB.M!=N || matB.N!=N){
608 std::cerr <<
"[ERROR] dgematrix::dggev"
609 <<
"(dgematrix&, vector<double>&, vector<double>&) " << std::endl
610 <<
"The matrix B is not a square matrix "
611 <<
"having the same size as \"this\" matrix." << std::endl
612 <<
"The B matrix is (" << matB.M <<
"x" << matB.N <<
")."
618 wr.resize(N); wi.resize(N);
619 char JOBVL(
'N'), JOBVR(
'N');
620 long LDA(N), LDB(N), LDVL(1), LDVR(1), LWORK(8*N), INFO(1);
621 double *BETA(
new double[N]), *VL(NULL), *VR(NULL),
622 *WORK(
new double[LWORK]);
623 dggev_(JOBVL, JOBVR, N,
Array, LDA, matB.Array, LDB, &wr[0], &wi[0],
624 BETA, VL, LDVL, VR, LDVR, WORK, LWORK, INFO);
625 delete [] WORK;
delete [] VL;
delete [] VR;
628 for(
long i=0;
i<N;
i++){ wr[
i]/=BETA[
i]; wi[
i]/=BETA[
i]; }
632 std::cerr <<
"[WARNING] dgematrix::dgeev"
633 <<
"(dgematrix&, vector<double>&, vector<double>&)"
635 <<
"Serious trouble happend. INFO = " << INFO <<
"."
650 std::vector<double>& wr, std::vector<double>& wi,
651 std::vector<dcovector>& vrr,
652 std::vector<dcovector>& vri)
655 std::cerr <<
"# [MARK] dgematrix::dggev(dgematrix&, std::vector<double>&, std::vector<double>&, std::vector<dcovector>&, std::vector<dcovector>&)"
661 std::cerr <<
"[ERROR] dgematrix::dggev"
662 <<
"(dgematrix&, vector<double>&, vector<double>&, "
663 <<
"vector<dcovector>&, vector<dcovector>&)" << std::endl
664 <<
"This matrix is not a square matrix." << std::endl
665 <<
"This matrix is (" << M <<
"x" << N <<
")." << std::endl;
668 if(matB.M!=N || matB.N!=N){
669 std::cerr <<
"[ERROR] dgematrix::dggev"
670 <<
"(dgematrix&, vector<double>&, vector<double>&, "
671 <<
"vector<dcovector>&, vector<dcovector>&)" << std::endl
672 <<
"The matrix B is not a square matrix "
673 <<
"having the same size as \"this\" matrix." << std::endl
674 <<
"The B matrix is (" << matB.M <<
"x" << matB.N <<
")."
680 wr.resize(N); wi.resize(N); vrr.resize(N); vri.resize(N);
681 for(
long i=0;
i<N;
i++){ vrr[
i].resize(N); vri[
i].resize(N); }
683 char JOBVL(
'N'), JOBVR(
'V');
684 long LDA(N), LDB(N), LDVL(1), LDVR(N), LWORK(8*N), INFO(1);
685 double *BETA(
new double[N]), *VL(NULL), *WORK(
new double[LWORK]);
686 dggev_(JOBVL, JOBVR, N,
Array, LDA, matB.Array, LDB, &wr[0], &wi[0],
687 BETA, VL, LDVL, VR.Array, LDVR, WORK, LWORK, INFO);
688 delete [] WORK;
delete [] VL;
691 for(
long i=0;
i<N;
i++){ wr[
i]/=BETA[
i]; wi[
i]/=BETA[
i]; }
695 for(
long j=0; j<N; j++){
696 if(fabs(wi[j])<1e-10){
697 for(
long i=0;
i<N;
i++){
698 vrr[j](
i) = VR(
i,j); vri[j](
i) = 0.0;
702 for(
long i=0;
i<N;
i++){
703 vrr[j](
i) = VR(
i,j); vri[j](
i) = VR(
i,j+1);
704 vrr[j+1](
i) = VR(
i,j); vri[j+1](
i) =-VR(
i,j+1);
711 std::cerr <<
"[WARNING] dgematrix::dggev"
712 <<
"(dgematrix&, vector<double>&, vector<double>&, "
713 <<
"vector<dcovector>&, vector<dcovector>&)" << std::endl
714 <<
"Serious trouble happend. INFO = " << INFO <<
"."
729 std::vector<double>& wr, std::vector<double>& wi,
730 std::vector<drovector>& vlr,
731 std::vector<drovector>& vli)
734 std::cerr <<
"# [MARK] dgematrix::dggev(dgematrix&, std::vector<double>&, std::vector<double>&, std::vector<drovector>&, std::vector<drovector>&)"
740 std::cerr <<
"[ERROR] dgematrix::dggev"
741 <<
"(dgematrix&, vector<double>&, vector<double>&, "
742 <<
"vector<drovector>&, vector<drovector>&)" << std::endl
743 <<
"This matrix is not a square matrix." << std::endl
744 <<
"This matrix is (" << M <<
"x" << N <<
")." << std::endl;
747 if(matB.M!=N || matB.N!=N){
748 std::cerr <<
"[ERROR] dgematrix::dggev"
749 <<
"(dgematrix&, vector<double>&, vector<double>&, "
750 <<
"vector<drovector>&, vector<drovector>&)" << std::endl
751 <<
"The matrix B is not a square matrix "
752 <<
"having the same size as \"this\" matrix." << std::endl
753 <<
"The B matrix is (" << matB.M <<
"x" << matB.N <<
")."
759 wr.resize(N); wi.resize(N); vlr.resize(N); vli.resize(N);
760 for(
long i=0;
i<N;
i++){ vlr[
i].resize(N); vli[
i].resize(N); }
762 char JOBVL(
'V'), JOBVR(
'N');
763 long LDA(N), LDB(N), LDVL(N), LDVR(1), LWORK(8*N), INFO(1);
764 double *BETA(
new double[N]), *VR(NULL), *WORK(
new double[LWORK]);
765 dggev_(JOBVL, JOBVR, N,
Array, LDA, matB.Array, LDB, &wr[0], &wi[0],
766 BETA, VL.Array, LDVL, VR, LDVR, WORK, LWORK, INFO);
767 delete [] WORK;
delete [] VR;
770 for(
long i=0;
i<N;
i++){ wr[
i]/=BETA[
i]; wi[
i]/=BETA[
i]; }
774 for(
long j=0; j<N; j++){
775 if(fabs(wi[j])<1e-10){
776 for(
long i=0;
i<N;
i++){
777 vlr[j](
i) = VL(
i,j); vli[j](
i) = 0.0;
781 for(
long i=0;
i<N;
i++){
782 vlr[j](
i) = VL(
i,j); vli[j](
i) =-VL(
i,j+1);
783 vlr[j+1](
i) = VL(
i,j); vli[j+1](
i) = VL(
i,j+1);
790 std::cerr <<
"[WARNING] dgematrix::dggev"
791 <<
"(dgematrix&, vector<double>&, vector<double>&, "
792 <<
"vector<drovector>&, vector<drovector>&)" << std::endl
793 <<
"Serious trouble happend. INFO = " << INFO <<
"."
815 std::cerr <<
"# [MARK] dgematrix::dgesvd(dcovector&, dgematrix&, dgematrix&)"
819 char JOBU(
'A'), JOBVT(
'A');
820 long LDA(M), LDU(M), LDVT(N),
821 LWORK(max(3*min(M,N)+max(M,N),5*min(M,N))), INFO(1);
822 double *WORK(
new double[LWORK]);
825 dgesvd_(JOBU, JOBVT, M, N,
Array, LDA, S.Array, U.Array,
826 LDU, VT.Array, LDVT, WORK, LWORK, INFO);
830 std::cerr <<
"[WARNING] dgematrix::dgesvd"
831 <<
"(dceovector&, dgematrix&, dcovector&) "
832 <<
"Serious trouble happend. INFO = " << INFO <<
"."
void clear()
Definition: dgematrix-misc.hpp:3
long dgelss(dcovector &, dcovector &, long &, const double)
Definition: dgematrix-lapack.hpp:304
long dgels(dgematrix &)
Definition: dgematrix-lapack.hpp:76
void resize(const long &)
Definition: dcovector-misc.hpp:93
void resize(const long &, const long &)
Definition: dgematrix-misc.hpp:126
friend void swap(dgematrix &, dgematrix &)
Definition: dgematrix-misc.hpp:154
double * Array
1D Array to store vector data
Definition: _drovector.hpp:8
Real Double-precision General Dence Matrix Class.
Definition: dgematrix.hpp:3
void zero()
Definition: drovector-misc.hpp:23
void clear()
Definition: dcovector-misc.hpp:3
Real Double-precision Row Vector Class.
Definition: drovector.hpp:3
long dgesv(dgematrix &)
Definition: dgematrix-lapack.hpp:6
Real Double-precision Column Vector Class.
Definition: dcovector.hpp:3
long dgeev(std::vector< double > &, std::vector< double > &)
Definition: dgematrix-lapack.hpp:410
friend _dgematrix i(const dgematrix &)
Definition: dgematrix-calc.hpp:21
void resize(const long &)
Definition: drovector-misc.hpp:93
long dgesvd(dcovector &, dgematrix &, dgematrix &)
Definition: dgematrix-lapack.hpp:812
long dggev(dgematrix &, std::vector< double > &, std::vector< double > &)
Definition: dgematrix-lapack.hpp:591