00001 #include "gmatrix.h"
00002 #include "precond.h"
00003
00004 gmatrix::gmatrix ()
00005 {
00006
00007 n=0;
00008
00009 zero=1.0e-30;
00010 fprintf (stdout,"\n computer zero in gmatrix is defined as %le",zero);
00011
00012
00013 nicg=0;
00014
00015 errcg=0.0;
00016
00017 anicg=0;
00018
00019 aerrcg=0.0;
00020
00021
00022 omega=0.0;
00023
00024 indegamma=0.0;
00025
00026 limit=0.0;
00027
00028
00029
00030 dm=NULL;
00031
00032 sky=NULL;
00033
00034 dsky=NULL;
00035
00036 cr=NULL;
00037
00038 scr=NULL;
00039
00040 em=NULL;
00041
00042
00043 sdirect = NULL;
00044
00045 auxsd=NULL;
00046
00047 bsize=0;
00048
00049 pmat=0;
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059 }
00060
00061 gmatrix::~gmatrix ()
00062 {
00063 delete dm;
00064 delete sky;
00065 delete dsky;
00066 delete cr;
00067 delete scr;
00068 delete em;
00069
00070 delete sdirect;
00071 delete auxsd;
00072
00073
00074
00075
00076
00077 }
00078
00079
00080
00081
00082
00083
00084 void gmatrix::alloc ()
00085 {
00086 switch (ts){
00087 case dense_matrix:{
00088 if (dm==NULL)
00089 dm = new densemat ();
00090 break;
00091 }
00092 case skyline_matrix:{
00093 if (sky==NULL)
00094 sky = new skyline ();
00095 break;
00096 }
00097 case double_skyline:{
00098 if (dsky==NULL)
00099 dsky = new dskyline ();
00100 break;
00101 }
00102 case compressed_rows:{
00103 if (cr==NULL)
00104 cr = new comprow ();
00105 break;
00106 }
00107 case symm_comp_rows:{
00108 if (scr==NULL)
00109 scr = new symcomprow ();
00110 break;
00111 }
00112 case element_matrices:{
00113 if (em==NULL)
00114 em = new elemmat ();
00115 break;
00116 }
00117
00118
00119 case lapack_stor:{
00120
00121 print_err("LAPACK solver is not supported at this time",__FILE__,__LINE__,__func__);
00122 abort ();
00123 break;
00124 }
00125
00126 case petsc:{
00127
00128 print_err("PETSC solver is not supported at this time",__FILE__,__LINE__,__func__);
00129 abort ();
00130 break;
00131 }
00132
00133 case spdirect_stor_cr:{
00134 if (cr==NULL)
00135 cr = new comprow ();
00136 if (sdirect==NULL)
00137 sdirect = new DSSolver();
00138
00139 break;
00140 }
00141
00142 case spdirect_stor_scr:{
00143 if (scr==NULL)
00144 scr = new symcomprow ();
00145 if (sdirect==NULL)
00146 sdirect = new DSSolver();
00147 break;
00148 }
00149
00150 default:{
00151 print_err("unknown type of matrix storage is required",__FILE__,__LINE__,__func__);
00152 }
00153 }
00154
00155
00156 pmat=0;
00157 }
00158
00159
00160
00161
00162 void gmatrix::dealloc ()
00163 {
00164 switch (ts){
00165 case dense_matrix:{
00166 delete dm;
00167 break;
00168 }
00169 case skyline_matrix:{
00170 delete sky;
00171 break;
00172 }
00173 case double_skyline:{
00174 delete dsky;
00175 break;
00176 }
00177 case compressed_rows:{
00178 delete cr;
00179 break;
00180 }
00181 case symm_comp_rows:{
00182 delete scr;
00183 break;
00184 }
00185 case element_matrices:{
00186 delete em;
00187 break;
00188 }
00189
00190
00191 case lapack_stor:{
00192
00193 print_err("LAPACK solver is not supported at this time",__FILE__,__LINE__,__func__);
00194 abort ();
00195 break;
00196 }
00197
00198 case petsc:{
00199
00200 print_err("PETSC solver is not supported at this time",__FILE__,__LINE__,__func__);
00201 abort ();
00202 break;
00203 }
00204
00205 case spdirect_stor_cr:{
00206 delete cr;
00207 delete sdirect;
00208
00209
00210 break;
00211 }
00212
00213 case spdirect_stor_scr:{
00214 delete scr;
00215 delete sdirect;
00216 break;
00217 }
00218
00219 default:{
00220 print_err("unknown type of matrix storage is required",__FILE__,__LINE__,__func__);
00221 }
00222 }
00223
00224 }
00225
00226
00227
00228
00229
00230
00231
00232
00233
00234
00235
00236 void gmatrix::initiate (gtopology *top,long ndof,storagetype ats,long mespr,FILE *out)
00237 {
00238
00239 n=ndof;
00240
00241 ts = ats;
00242
00243 alloc ();
00244
00245 switch (ts){
00246 case dense_matrix:{
00247 dm->initiate (n,mespr);
00248 break;
00249 }
00250 case skyline_matrix:{
00251 sky->initiate (top,n,mespr);
00252 break;
00253 }
00254 case double_skyline:{
00255 dsky->initiate (top,n,mespr);
00256 break;
00257 }
00258 case compressed_rows:{
00259 cr->initiate (top,n,mespr);
00260 break;
00261 }
00262 case symm_comp_rows:{
00263 scr->initiate (top,n,mespr);
00264 break;
00265 }
00266 case element_matrices:{
00267 em->initiate (top,n,mespr);
00268 break;
00269 }
00270
00271 case lapack_stor:{
00272
00273
00274
00275 print_err("LAPACK solver is not supported at this time",__FILE__,__LINE__,__func__);
00276 abort ();
00277 break;
00278 }
00279
00280 case petsc:{
00281
00282
00283
00284 print_err("PETSC solver is not supported at this time",__FILE__,__LINE__,__func__);
00285 abort ();
00286 break;
00287 }
00288
00289 case spdirect_stor_cr:{
00290 cr->initiate (top,n,mespr);
00291 break;
00292 }
00293 case spdirect_stor_scr:{
00294 scr->initiate (top,n,mespr);
00295 break;
00296 }
00297
00298 default:{
00299 print_err("unknown type of matrix storage is required",__FILE__,__LINE__,__func__);
00300 }
00301 }
00302
00303 if (ts == spdirect_stor_cr || ts == spdirect_stor_scr){
00304 if (sdirect->IsInitialized () == 0){
00305
00306
00307
00308
00309
00310
00311
00312
00313 switch (tlinsol){
00314 case spdirldl:{
00315 sdirect->Initialize(1,eDSSFactorizationLDLT);
00316 break;
00317 }
00318 case spdirlu:{
00319 sdirect->Initialize(1,eDSSFactorizationLU);
00320 break;
00321 }
00322 case spdirll:{
00323 sdirect->Initialize(1,eDSSFactorizationLLT);
00324 break;
00325 }
00326 case cg:{
00327 if (tprec == sparseindec)
00328 sdirect->Initialize(1,eDSSFactorizationLDLTIncomplete);
00329 break;
00330 }
00331 case conden:{
00332 switch (tsol){
00333 case spdirldl:{
00334 sdirect->Initialize(1,eDSSFactorizationLDLT);
00335 break;
00336 }
00337 case spdirlu:{
00338 sdirect->Initialize(1,eDSSFactorizationLU);
00339 break;
00340 }
00341 case spdirll:{
00342 sdirect->Initialize(1,eDSSFactorizationLLT);
00343 break;
00344 }
00345 case cg:{
00346 if (tprec == sparseindec)
00347 sdirect->Initialize(1,eDSSFactorizationLDLTIncomplete);
00348 break;
00349 }
00350 default:{
00351 print_err("unknown type of sparse solver is required",__FILE__,__LINE__,__func__);
00352 }
00353 }
00354
00355 break;
00356 }
00357
00358 default:{
00359 print_err("unknown type of sparse solver is required",__FILE__,__LINE__,__func__);
00360 }
00361 }
00362
00363
00364
00365 if (top->gnodes == NULL){
00366 bsize = (unsigned char) 6;
00367 }
00368 else{
00369 bsize = (unsigned char) top->gnodes[0].ndofn;
00370 }
00371 if (bsize<1){
00372 print_err("wrong block size",__FILE__,__LINE__,__func__);
00373 abort ();
00374 }
00375
00376
00377
00378
00379
00380
00381
00382
00383 if (ts == spdirect_stor_cr){
00384 SparseMatrixF sm((unsigned long)n,cr->a,(unsigned long*)cr->ci,(unsigned long*)cr->adr,0,0,1,0);
00385
00386 sdirect->SetMatrixPattern (&sm,bsize);
00387 }
00388 else{
00389 SparseMatrixF sm((unsigned long)n,scr->a,(unsigned long*)scr->ci,(unsigned long*)scr->adr,0,0,1,1);
00390
00391 sdirect->SetMatrixPattern (&sm,bsize);
00392 }
00393
00394 if (tlinsol==conden){
00395
00396 auxdatsparsesolver (top,out);
00397
00398 if (top->nn != 0){
00399 sdirect->LoadMCN (top->nn,bsize,auxsd,1);
00400 }
00401
00402 }
00403
00404
00405 }
00406 else{
00407
00408 sdirect->LoadZeros ();
00409 }
00410 }
00411
00412 }
00413
00414
00415
00416
00417
00418
00419
00420
00421
00422 void gmatrix::setval (slesolv *ssle)
00423 {
00424
00425 tlinsol=ssle->tlinsol;
00426
00427 tsol=ssle->tsol;
00428
00429 tprec=ssle->prec.pt;
00430
00431
00432 nicg=ssle->ni;
00433
00434 errcg=ssle->err;
00435
00436 iv=ssle->iv;
00437
00438
00439
00440 }
00441
00442
00443
00444
00445
00446
00447
00448
00449
00450
00451
00452
00453 void gmatrix::localize (matrix &lm,ivector &cn,long eid)
00454 {
00455 if (lm.m != cn.n){
00456
00457 }
00458
00459 switch (ts){
00460 case dense_matrix:{
00461 dm->localize (lm,cn.a);
00462 break;
00463 }
00464 case skyline_matrix:{
00465 sky->localize (lm,cn.a);
00466 break;
00467 }
00468 case double_skyline:{
00469 dsky->localize (lm,cn.a);
00470 break;
00471 }
00472 case compressed_rows:{
00473 cr->localize (lm,cn.a);
00474 break;
00475 }
00476 case symm_comp_rows:{
00477 scr->localize (lm,cn.a);
00478 break;
00479 }
00480 case element_matrices:{
00481 em->localize (lm,cn.a,eid);
00482 break;
00483 }
00484
00485 case lapack_stor:{
00486
00487 print_err("LAPACK solver is not supported at this time",__FILE__,__LINE__,__func__);
00488 abort ();
00489 break;
00490 }
00491 case petsc:{
00492
00493 print_err("PETSC solver is not supported at this time",__FILE__,__LINE__,__func__);
00494 abort ();
00495 break;
00496 }
00497
00498 case spdirect_stor_cr:{
00499 cr->localize (lm,cn.a);
00500 break;
00501 }
00502 case spdirect_stor_scr:{
00503 scr->localize (lm,cn.a);
00504 break;
00505 }
00506 default:{
00507 print_err("unknown type of matrix storage is required",__FILE__,__LINE__,__func__);
00508 }
00509 }
00510 }
00511
00512
00513
00514
00515
00516
00517
00518
00519
00520
00521
00522
00523 void gmatrix::localized (double *lm,long *cn,long nc,long eid)
00524 {
00525 switch (ts){
00526 case dense_matrix:{
00527 dm->localized (lm,cn,nc);
00528 break;
00529 }
00530 case skyline_matrix:{
00531 sky->localized (lm,cn,nc);
00532 break;
00533 }
00534 case double_skyline:{
00535 dsky->localized (lm,cn,nc);
00536 break;
00537 }
00538 case compressed_rows:{
00539 cr->localized (lm,cn,nc);
00540 break;
00541 }
00542 case symm_comp_rows:{
00543 scr->localized (lm,cn,nc);
00544 break;
00545 }
00546 case element_matrices:{
00547 em->localized (lm,cn,nc,eid);
00548 break;
00549 }
00550 case lapack_stor:{
00551
00552 print_err("LAPACK solver is not supported at this time",__FILE__,__LINE__,__func__);
00553 abort ();
00554 break;
00555 }
00556 case petsc:{
00557
00558 print_err("PETSC solver is not supported at this time",__FILE__,__LINE__,__func__);
00559 abort ();
00560 break;
00561 }
00562 case spdirect_stor_cr:{
00563 cr->localized (lm,cn,nc);
00564 break;
00565 }
00566 case spdirect_stor_scr:{
00567 scr->localized (lm,cn,nc);
00568 break;
00569 }
00570 default:{
00571 print_err("unknown type of matrix storage is required",__FILE__,__LINE__,__func__);
00572 }
00573 }
00574 }
00575
00576
00577
00578
00579
00580
00581
00582
00583
00584
00585
00586
00587 void gmatrix::glocalize (matrix &lm,ivector &rcn,ivector &ccn)
00588 {
00589 switch (ts){
00590 case dense_matrix:{
00591 dm->glocalize (lm,rcn.a,ccn.a);
00592 break;
00593 }
00594 case skyline_matrix:{
00595 sky->glocalize (lm,rcn.a,ccn.a);
00596 break;
00597 }
00598 case double_skyline:{
00599
00600 break;
00601 }
00602 case compressed_rows:{
00603
00604 break;
00605 }
00606 case symm_comp_rows:{
00607
00608 break;
00609 }
00610 default:{
00611 print_err("unknown type of matrix storage is required",__FILE__,__LINE__,__func__);
00612 }
00613 }
00614 }
00615
00616
00617
00618
00619
00620
00621
00622
00623
00624
00625
00626 void gmatrix::mult_localize (long nm,long *ncn1,long *ncn2,long *mcn)
00627 {
00628 switch (ts){
00629 case dense_matrix:{
00630 dm->mult_localize (nm,ncn1,ncn2,mcn);
00631 break;
00632 }
00633 case skyline_matrix:{
00634 sky->mult_localize (nm,ncn1,ncn2,mcn);
00635 break;
00636 }
00637 case double_skyline:{
00638 dsky->mult_localize (nm,ncn1,ncn2,mcn);
00639 break;
00640 }
00641 default:{
00642 print_err("wrong storage type of matrix is required",__FILE__,__LINE__,__func__);
00643 }
00644 }
00645 }
00646
00647
00648
00649
00650
00651
00652
00653
00654
00655
00656 void gmatrix::prepmat (double limit,long mespr)
00657 {
00658 switch (ts){
00659 case dense_matrix:{
00660 dm->decompid=0;
00661 break;
00662 }
00663 case skyline_matrix:{
00664 sky->setnotfact ();
00665 break;
00666 }
00667 case double_skyline:{
00668 dsky->setnotfact ();
00669 break;
00670 }
00671 case compressed_rows:{
00672 cr->decompid=0;
00673 long rejected=0;
00674
00675 if (mespr==1) fprintf (stdout,"\n number of rejected entries %ld",rejected);
00676 if (mespr==1) fprintf (stdout,"\n number of matrix entries %ld",cr->negm);
00677 break;
00678 }
00679 case symm_comp_rows:{
00680 scr->decompid=0;
00681 long rejected=0;
00682
00683 if (mespr==1) fprintf (stdout,"\n number of rejected entries %ld",rejected);
00684 if (mespr==1) fprintf (stdout,"\n number of matrix entries %ld",scr->negm);
00685 break;
00686 }
00687
00688 case lapack_stor:{
00689
00690
00691
00692
00693
00694
00695
00696
00697 print_err("LAPACK solver is not supported at this time",__FILE__,__LINE__,__func__);
00698 abort ();
00699 break;
00700 }
00701
00702 case petsc:{
00703
00704
00705
00706
00707
00708
00709
00710
00711
00712
00713
00714
00715
00716
00717
00718
00719
00720
00721
00722
00723
00724
00725
00726
00727
00728
00729
00730
00731
00732 print_err("PETSC solver is not supported at this time",__FILE__,__LINE__,__func__);
00733 abort ();
00734 break;
00735 }
00736
00737
00738 case spdirect_stor_cr:{
00739
00740 break;
00741 }
00742 case spdirect_stor_scr:{
00743
00744 break;
00745 }
00746
00747
00748 case element_matrices:{
00749 em->decompid=0;
00750 break;
00751 }
00752
00753 default:{
00754 print_err("unknown type of matrix storage is required",__FILE__,__LINE__,__func__);
00755 }
00756 }
00757
00758 }
00759
00760
00761
00762
00763
00764
00765
00766
00767
00768
00769 void gmatrix::prepmat2 (gtopology *top,FILE *out)
00770 {
00771 switch (ts){
00772 case dense_matrix:{
00773 break;
00774 }
00775 case skyline_matrix:{
00776 break;
00777 }
00778 case double_skyline:{
00779 break;
00780 }
00781 case compressed_rows:{
00782 break;
00783 }
00784 case symm_comp_rows:{
00785 break;
00786 }
00787
00788 case lapack_stor:{
00789 print_err("LAPACK solver is not supported at this time",__FILE__,__LINE__,__func__);
00790 abort ();
00791 break;
00792 }
00793 case petsc:{
00794 print_err("PETSC solver is not supported at this time",__FILE__,__LINE__,__func__);
00795 abort ();
00796 break;
00797 }
00798
00799
00800 case spdirect_stor_cr:{
00801 if (sdirect->IsAllocated()==0){
00802
00803
00804
00805
00806 sdirect->PreFactorize ();
00807 }
00808
00809
00810
00811
00812
00813
00814
00815
00816 SparseMatrixF sm((unsigned long)n,cr->a,(unsigned long*)cr->ci,(unsigned long*)cr->adr,0,0,1,0);
00817 sdirect->LoadNumbers (&sm);
00818
00819
00820
00821 break;
00822 }
00823
00824 case spdirect_stor_scr:{
00825 if (sdirect->IsAllocated()==0){
00826
00827
00828
00829
00830 sdirect->PreFactorize ();
00831 }
00832
00833
00834
00835
00836
00837
00838
00839
00840 SparseMatrixF sm((unsigned long)n,scr->a,(unsigned long*)scr->ci,(unsigned long*)scr->adr,0,0,1,1);
00841 sdirect->LoadNumbers (&sm);
00842
00843
00844
00845 break;
00846 }
00847
00848
00849 case element_matrices:{
00850 break;
00851 }
00852
00853 default:{
00854 print_err("unknown type of matrix storage is required",__FILE__,__LINE__,__func__);
00855 }
00856 }
00857
00858
00859 pmat=1;
00860 }
00861
00862
00863
00864
00865
00866
00867
00868
00869
00870
00871
00872 void gmatrix::auxdatsparsesolver (gtopology *top,FILE *out)
00873 {
00874 long i,j,k,m,nn,ndofn,rndofn;
00875
00876
00877 fprintf (out,"\n\n\n\n kontrola Richarda \n\n");
00878 for (i=0;i<top->nn;i++){
00879 fprintf (out,"\n %ld %ld",i,top->gnodes[i].ai);
00880 }
00881
00882
00883
00884
00885
00886
00887
00888
00889
00890 if (top->gnodes != NULL){
00891
00892 nn = top->nn;
00893 ndofn = top->gnodes[0].ndofn;
00894
00895
00896
00897
00898 auxsd = new long [nn*ndofn];
00899
00900 k=0;
00901 for (i=0;i<nn;i++){
00902 rndofn=top->gnodes[i].ndofn;
00903 if (rndofn!=ndofn){
00904 print_err("wrong number of degrees of freedom of node is required",__FILE__,__LINE__,__func__);
00905 abort ();
00906 }
00907
00908 if (top->gnodes[i].ai>-1){
00909
00910 for (j=0;j<ndofn;j++){
00911 m=top->give_dof (i,j);
00912 if (m>-1){
00913
00914 auxsd[k] = -1 - top->give_dof (i,j);
00915 k++;
00916 }
00917 if (m<0){
00918
00919 auxsd[k] = -1;
00920 k++;
00921 }
00922 }
00923 }
00924 else{
00925
00926 for (j=0;j<ndofn;j++){
00927 m=top->give_dof (i,j);
00928 if (m>-1){
00929
00930 auxsd[k] = m-1;
00931 k++;
00932 }
00933 if (m<0){
00934
00935 auxsd[k] = -1;
00936 k++;
00937 }
00938 }
00939 }
00940 }
00941
00942 }
00943
00944 fprintf (out,"\n\n\n\n kontrola Richarda \n\n");
00945 k=0;
00946 for (i=0;i<top->nn;i++){
00947 fprintf (out,"\n uzel %5ld",i);
00948 for (j=0;j<ndofn;j++){
00949 fprintf (out," %ld",auxsd[k]);
00950 k++;
00951 }
00952 }
00953 fprintf (out,"\n\n\n");
00954
00955 }
00956
00957
00958
00959
00960
00961
00962
00963
00964
00965
00966
00967
00968 void gmatrix::solve_system (gtopology *top,precond &prec,double *lhs,double *rhs,FILE *out)
00969 {
00970 time_t bt,et;
00971 bt = time (NULL);
00972
00973 if (pmat==0)
00974 prepmat2 (top,out);
00975
00976 switch (tlinsol){
00977
00978 case gauss_elim:{
00979 switch (ts){
00980 case dense_matrix:{
00981 dm->gemp (lhs,rhs,1,zero,2);
00982 break;
00983 }
00984 default:{
00985 print_err("wrong storage type of matrix is required",__FILE__,__LINE__,__func__);
00986 }
00987 }
00988 break;
00989 }
00990
00991 case ll:{
00992 switch (ts){
00993 case dense_matrix:{
00994 dm->ll (lhs,rhs,zero,1);
00995 break;
00996 }
00997 default:{
00998 print_err("wrong storage type of matrix is required",__FILE__,__LINE__,__func__);
00999 }
01000 }
01001 break;
01002 }
01003
01004 case ill:{
01005 switch (ts){
01006 case dense_matrix:{
01007 dm->ill (lhs,rhs,zero,zero,1);
01008 break;
01009 }
01010 default:{
01011 print_err("wrong storage type of matrix is required",__FILE__,__LINE__,__func__);
01012 }
01013 }
01014 break;
01015 }
01016
01017 case ldl:{
01018 switch (ts){
01019 case dense_matrix:{
01020 if (dm->decomp()==0){
01021 dm->ll (NULL,NULL,zero,2);
01022 }
01023 dm->ll (lhs,rhs,zero,3);
01024 break;
01025 }
01026 case skyline_matrix:{
01027 if (sky->decomp()==0){
01028 sky->ldl_sky (NULL,NULL,zero,2);
01029 }
01030 sky->ldl_sky (lhs,rhs,zero,3);
01031 break;
01032 }
01033 default:{
01034 print_err("wrong storage type of matrix is required",__FILE__,__LINE__,__func__);
01035 }
01036 }
01037 break;
01038 }
01039
01040 case lu:{
01041 switch (ts){
01042 case dense_matrix:{
01043 if (dm->decomp()==0){
01044 dm->lu (lhs,rhs,zero,2);
01045 }
01046 dm->lu (lhs,rhs,zero,3);
01047 break;
01048 }
01049 case double_skyline:{
01050 if (dsky->decomp()==0){
01051 dsky->lu_dsky (NULL,NULL,zero,2);
01052 }
01053 dsky->lu_dsky (lhs,rhs,zero,3);
01054 break;
01055 }
01056 default:{
01057 print_err("wrong storage type of matrix is required",__FILE__,__LINE__,__func__);
01058 }
01059 }
01060 break;
01061 }
01062
01063 case conden:{
01064 double *condvect,*av;
01065 densemat adm;
01066
01067 condvect=new double[nbdof];
01068 av=new double[nbdof];
01069
01070 adm.alloc(nbdof);
01071
01072 condense (top,adm.a,condvect,lhs,rhs,nbdof,1,out);
01073
01074 adm.gemp (av,condvect,1,1.0e-10,1);
01075
01076 copyv (av,condvect,nbdof);
01077
01078 condense (top,adm.a,condvect,lhs,rhs,nbdof,2,out);
01079
01080 delete [] condvect;
01081 delete [] av;
01082 adm.dealloc ();
01083
01084 break;
01085 }
01086
01087 case cg:{
01088 switch (ts){
01089 case dense_matrix:{
01090 switch (tprec){
01091 case noprecond:{
01092 dm->cg (lhs,rhs,nicg,errcg,anicg,aerrcg,zero,iv);
01093 break;
01094 }
01095 case incomdec:{
01096 fprintf (stdout,"\n incomplete decomposition of matrix");
01097 densemat dmid;
01098 dmid.copy (dm);
01099
01100 if (dmid.decomp()==0){
01101 dmid.ill (lhs,rhs,zero,zero,2);
01102
01103 }
01104
01105
01106 dm->cg_prec_new (prec,lhs,rhs,nicg,errcg,anicg,aerrcg,zero,iv);
01107 break;
01108 }
01109 default:{
01110 print_err("unknown preconditioner is required",__FILE__,__LINE__,__func__);
01111 }
01112 }
01113
01114 break;
01115 }
01116 case compressed_rows:{
01117 cr->cg_prec (prec,lhs,rhs,nicg,errcg,anicg,aerrcg,zero,iv);
01118
01119
01120
01121 break;
01122 }
01123 case symm_comp_rows:{
01124 switch (tprec){
01125 case noprecond:{
01126 scr->cg (lhs,rhs,nicg,errcg,anicg,aerrcg,zero,iv);
01127 break;
01128 }
01129 case diagprec:{
01130 scr->cg_prec (lhs,rhs,nicg,errcg,anicg,aerrcg,zero,iv,tprec,0.0,sdirect);
01131 break;
01132 }
01133 case ssorprec:{
01134 scr->cg_prec (lhs,rhs,nicg,errcg,anicg,aerrcg,zero,iv,tprec,omega,sdirect);
01135 break;
01136 }
01137 case incomdec:{
01138 fprintf (stdout,"\n incomplete decomposition of matrix");
01139 if (scr->decomp()==0){
01140 scr->incomplete_ldl (indegamma);
01141 scr->changedecomp ();
01142 }
01143
01144 scr->cg_prec (lhs,rhs,nicg,errcg,anicg,aerrcg,zero,iv,tprec,0.0,sdirect);
01145 break;
01146 }
01147 case sparseindec:{
01148 fprintf (stdout,"\n incomplete sparse decomposition of matrix");
01149
01150 if (sdirect->IsFactorized()==0){
01151 sdirect->ReFactorize ();
01152 }
01153
01154 scr->cg_prec (lhs,rhs,nicg,errcg,anicg,aerrcg,zero,iv,tprec,0.0,sdirect);
01155 break;
01156 }
01157
01158 default:{
01159 fprintf (stderr,"\n\n unknown preconditioner type is required in function cgsol (file %s, line %d).\n",__FILE__,__LINE__);
01160 }
01161 }
01162
01163 break;
01164 }
01165 case element_matrices:{
01166 em->cg (lhs,rhs,nicg,errcg,anicg,aerrcg,zero,iv);
01167 break;
01168 }
01169
01170 default:{
01171 print_err("wrong storage type of matrix is required",__FILE__,__LINE__,__func__);
01172 }
01173 }
01174 break;
01175 }
01176
01177 case bicg:{
01178 switch (ts){
01179 case double_skyline:{
01180 dsky->bicg (lhs,rhs,nicg,errcg,anicg,aerrcg,zero,iv);
01181 break;
01182 }
01183 case compressed_rows:{
01184 cr->bicg (lhs,rhs,nicg,errcg,anicg,aerrcg,zero,iv);
01185
01186 break;
01187 }
01188 default:{
01189 print_err("wrong storage type of matrix is required",__FILE__,__LINE__,__func__);
01190 }
01191 }
01192 break;
01193 }
01194
01195 case lapack_sol:{
01196
01197
01198
01199 print_err("LAPACK solver is not supported at this time",__FILE__,__LINE__,__func__);
01200 abort ();
01201 break;
01202 }
01203
01204 case spdirldl:
01205 case spdirlu:
01206 case spdirll:{
01207 if (sdirect->IsFactorized()==0){
01208 sdirect->ReFactorize();
01209 }
01210 sdirect->Solve(lhs,rhs);
01211 break;
01212 }
01213
01214
01215 case pardisolver:{
01216
01217
01218
01219
01220
01221
01222
01223
01224
01225
01226
01227
01228
01229
01230
01231 print_err("PARDISO solver is not supported at this time",__FILE__,__LINE__,__func__);
01232
01233 break;
01234 }
01235
01236 default:{
01237 print_err("unknown type of linear solveris required",__FILE__,__LINE__,__func__);
01238 }
01239 }
01240
01241 et = time (NULL);
01242
01243 fflush(stdout);
01244 }
01245
01246
01247
01248
01249
01250
01251 void gmatrix::decompose_matrix ()
01252 {
01253 switch (tlinsol){
01254
01255 case ldl:{
01256 switch (ts){
01257 case skyline_matrix:{
01258 if (sky->decomp()==0){
01259 sky->ldl_sky (NULL,NULL,zero,2);
01260 }
01261 break;
01262 }
01263 default:{
01264 print_err("wrong storage type of matrix is required",__FILE__,__LINE__,__func__);
01265 }
01266 }
01267 break;
01268 }
01269
01270 case lu:{
01271 switch (ts){
01272 case dense_matrix:{
01273 if (dm->decomp()==0){
01274 dm->lu (NULL,NULL,zero,2);
01275 }
01276 break;
01277 }
01278 case double_skyline:{
01279 if (dsky->decomp()==0){
01280 dsky->lu_dsky (NULL,NULL,zero,2);
01281 }
01282 break;
01283 }
01284 default:{
01285 print_err("wrong storage type of matrix is required",__FILE__,__LINE__,__func__);
01286 }
01287 }
01288 break;
01289 }
01290
01291 case lapack_sol:{
01292
01293
01294 print_err("LAPACK solver is not supported at this time",__FILE__,__LINE__,__func__);
01295 abort ();
01296 break;
01297 }
01298
01299 case spdirldl:
01300 case spdirlu:
01301 case spdirll:{
01302 if (sdirect->IsFactorized()==0){
01303 sdirect->ReFactorize();
01304 }
01305 break;
01306 }
01307
01308
01309 default:{
01310 print_err("unknown type of linear solver is required",__FILE__,__LINE__,__func__);
01311 }
01312 }
01313 }
01314
01315
01316
01317
01318
01319
01320
01321
01322
01323 void gmatrix::back_substitution (double *lhs,double *rhs)
01324 {
01325 switch (tlinsol){
01326
01327 case ldl:{
01328 switch (ts){
01329 case skyline_matrix:{
01330 if (sky->decomp()==0){
01331 sky->ldl_sky (NULL,NULL,zero,2);
01332 }
01333 sky->ldl_sky (lhs,rhs,zero,3);
01334 break;
01335 }
01336 default:{
01337 print_err("wrong storage type of matrix is required",__FILE__,__LINE__,__func__);
01338 }
01339 }
01340 break;
01341 }
01342
01343 case lu:{
01344 switch (ts){
01345 case dense_matrix:{
01346 if (dm->decomp()==0){
01347 dm->lu (lhs,rhs,zero,2);
01348 }
01349 dm->lu (lhs,rhs,zero,3);
01350 break;
01351 }
01352 case double_skyline:{
01353 if (dsky->decomp()==0){
01354 dsky->lu_dsky (NULL,NULL,zero,2);
01355 }
01356 dsky->lu_dsky (lhs,rhs,zero,3);
01357 break;
01358 }
01359 default:{
01360 print_err("wrong storage type of matrix is required",__FILE__,__LINE__,__func__);
01361 }
01362 }
01363 break;
01364 }
01365
01366 case lapack_sol:{
01367
01368 print_err("LAPACK solver is not supported at this time",__FILE__,__LINE__,__func__);
01369 abort ();
01370 break;
01371 }
01372
01373 case spdirldl:
01374 case spdirlu:
01375 case spdirll:{
01376 if (sdirect->IsFactorized()==0){
01377 sdirect->ReFactorize();
01378 }
01379 sdirect->Solve(lhs,rhs);
01380 break;
01381 }
01382 case pardisolver:{
01383
01384
01385
01386 print_err("PARDISO solver is not supported at this time",__FILE__,__LINE__,__func__);
01387
01388 break;
01389 }
01390 default:{
01391 print_err("unknown type of linear solver is required",__FILE__,__LINE__,__func__);
01392 }
01393 }
01394 }
01395
01396
01397
01398
01399
01400
01401
01402
01403 void gmatrix::incomplete_fact (double incompltresh)
01404 {
01405 double *x,*y;
01406 x=NULL;
01407 y=NULL;
01408
01409 switch (ts){
01410 case dense_matrix:{
01411 dm->ill (x,y,zero,incompltresh,2);
01412 break;
01413 }
01414 default:{
01415 print_err("unknown type of matrix storage is required",__FILE__,__LINE__,__func__);
01416 }
01417 }
01418
01419 }
01420
01421
01422
01423
01424
01425
01426
01427
01428
01429 void gmatrix::back_incomplete_fact (double *x,double *y,double incompltresh)
01430 {
01431
01432 switch (ts){
01433 case dense_matrix:{
01434 dm->ill (x,y,zero,incompltresh,3);
01435 break;
01436 }
01437 default:{
01438 print_err("unknown type of matrix storage is required",__FILE__,__LINE__,__func__);
01439 }
01440 }
01441
01442 }
01443
01444
01445
01446
01447
01448
01449
01450
01451
01452
01453
01454
01455
01456
01457
01458
01459
01460 void gmatrix::condense (gtopology *top,double *condmat,double *condvect,double *lhs,double *rhs,long nrdof,long tc,FILE *out)
01461 {
01462 long i,j;
01463
01464 if (pmat==0)
01465 prepmat2 (top,out);
01466
01467 switch (ts){
01468 case dense_matrix:{
01469 dm->gempkon (condmat,condvect,lhs,rhs,nrdof,zero,tc);
01470 break;
01471 }
01472 case skyline_matrix:{
01473 sky->ldlkon_sky (condmat,condvect,lhs,rhs,nrdof,tc);
01474 break;
01475 }
01476 case double_skyline:{
01477 dsky->lukon_dsky (condmat,condvect,lhs,rhs,zero,nrdof,tc);
01478 break;
01479 }
01480 case spdirect_stor_cr:
01481 case spdirect_stor_scr:{
01482
01483
01484 if (tc==2){
01485 j=0;
01486 for (i=top->nidof;i<n;i++){
01487 lhs[i]=condvect[j];
01488 j++;
01489 }
01490 }
01491
01492 sdirect->condense (condmat,lhs,rhs,tc);
01493
01494 if (tc==1){
01495 j=0;
01496 for (i=top->nidof;i<n;i++){
01497 condvect[j]=rhs[i];
01498 j++;
01499 }
01500 }
01501
01502 break;
01503 }
01504 default:{
01505 print_err("wrong storage type of matrix is required",__FILE__,__LINE__,__func__);
01506 }
01507 }
01508 }
01509
01510
01511
01512
01513
01514
01515
01516
01517
01518
01519
01520
01521
01522 void gmatrix::kernel (double *rbm,long &nse,long *se,long ense,double limit,long tc)
01523 {
01524 switch (ts){
01525 case dense_matrix:{
01526 dm->ker (rbm,nse,se,ense,limit);
01527 break;
01528 }
01529 case skyline_matrix:{
01530 sky->ker (rbm,nse,se,ense,limit,tc);
01531 break;
01532 }
01533 default:{
01534 print_err("wrong storage type of matrix is required",__FILE__,__LINE__,__func__);
01535 }
01536 }
01537 }
01538
01539
01540
01541
01542
01543
01544
01545 void gmatrix::ldl_feti (double *lhs,double *rhs,long nse,long *se,double zero)
01546 {
01547 switch (ts){
01548 case skyline_matrix:{
01549 sky->ldl_feti_sky (lhs,rhs,nse,se,zero);
01550 break;
01551 }
01552 default:{
01553 print_err("wrong storage type of matrix is required",__FILE__,__LINE__,__func__);
01554 }
01555 }
01556 }
01557
01558
01559
01560
01561
01562
01563
01564
01565 void gmatrix::gmxv (double *a,double *b)
01566 {
01567 switch (ts){
01568 case dense_matrix:{
01569 dm -> mxv_dm(a,b);
01570 break;
01571 }
01572 case skyline_matrix:{
01573 sky -> mxv_sky(a,b);
01574 break;
01575 }
01576 case double_skyline:{
01577 dsky -> mxv_dsky(a,b);
01578 break;
01579 }
01580 case compressed_rows:{
01581 cr -> mxv_cr(a,b);
01582 break;
01583 }
01584 case symm_comp_rows:{
01585 scr -> mxv_scr(a,b);
01586 break;
01587 }
01588 case element_matrices:{
01589 em -> mxv_em(a,b);
01590 break;
01591 }
01592 case spdirect_stor_scr:{
01593 scr -> mxv_scr(a,b);
01594 break;
01595 }
01596 case spdirect_stor_cr:{
01597 cr -> mxv_cr(a,b);
01598 break;
01599 }
01600 default:{
01601 print_err("unknown type of matrix storage is required",__FILE__,__LINE__,__func__);
01602 }
01603 }
01604
01605 }
01606
01607
01608
01609
01610
01611
01612
01613
01614
01615 void gmatrix::decompgmxv (double *a,double *b)
01616 {
01617 switch (ts){
01618 case skyline_matrix:{
01619 sky -> ldlmxv_sky (a,b);
01620 break;
01621 }
01622 case spdirect_stor_cr:{
01623 sdirect -> MulMatrixByVector (a,b);
01624 break;
01625 }
01626 case spdirect_stor_scr:{
01627 sdirect -> MulMatrixByVector (a,b);
01628 break;
01629 }
01630 default:{
01631 print_err("unknown type of matrix storage is required",__FILE__,__LINE__,__func__);
01632 }
01633 }
01634 }
01635
01636
01637
01638
01639
01640 void gmatrix::addgm (double a,gmatrix &gm)
01641 {
01642 switch (ts){
01643 case dense_matrix:{
01644 dm -> addmat_dm (a,*gm.dm);
01645 break;
01646 }
01647 case skyline_matrix:{
01648 sky -> addmat_sky (a,*gm.sky);
01649 break;
01650 }
01651 case double_skyline:{
01652 dsky -> addmat_dsky (a,*gm.dsky);
01653 break;
01654 }
01655 case compressed_rows:{
01656 cr -> addmat_cr (a,*gm.cr);
01657 break;
01658 }
01659 case symm_comp_rows:{
01660 scr -> addmat_scr (a,*gm.scr);
01661 break;
01662 }
01663 case element_matrices:{
01664 em -> addmat_em (a,*gm.em);
01665 break;
01666 }
01667 case spdirect_stor_scr:{
01668 scr -> addmat_scr (a,*gm.scr);
01669
01670 break;
01671 }
01672 case spdirect_stor_cr:{
01673 cr -> addmat_cr (a,*gm.cr);
01674 break;
01675 }
01676 default:{
01677 print_err("unknown type of matrix storage is required",__FILE__,__LINE__,__func__);
01678 }
01679 }
01680 }
01681
01682
01683
01684
01685
01686
01687
01688
01689 void gmatrix::scalgm (double a)
01690 {
01691 switch (ts){
01692 case dense_matrix:{
01693 dm -> scalmat_dm (a);
01694 break;
01695 }
01696 case skyline_matrix:{
01697 sky -> scalmat_sky (a);
01698 break;
01699 }
01700 case double_skyline:{
01701 dsky -> scalmat_dsky (a);
01702 break;
01703 }
01704 case compressed_rows:{
01705 cr -> scalmat_cr (a);
01706 break;
01707 }
01708 case symm_comp_rows:{
01709 scr -> scalmat_scr (a);
01710 break;
01711 }
01712 case element_matrices:{
01713 em -> scalmat_em (a);
01714 break;
01715 }
01716 case spdirect_stor_scr:{
01717 scr -> scalmat_scr (a);
01718
01719 break;
01720 }
01721 case spdirect_stor_cr:{
01722 cr -> scalmat_cr (a);
01723 break;
01724 }
01725 default:{
01726 print_err("unknown type of matrix storage is required",__FILE__,__LINE__,__func__);
01727 }
01728 }
01729 }
01730
01731
01732
01733
01734
01735
01736
01737
01738
01739 void gmatrix::copygm (gmatrix &gm)
01740 {
01741
01742 gm.ts=ts;
01743
01744 switch (ts){
01745 case dense_matrix:{
01746 if (gm.dm==NULL){
01747 gm.dm = new densemat ();
01748 }
01749 dm -> copy_dm (*gm.dm);
01750 break;
01751 }
01752 case skyline_matrix:{
01753 if (gm.sky==NULL){
01754 gm.sky = new skyline ();
01755 }
01756 sky -> copy_sky (*gm.sky);
01757 break;
01758 }
01759 case double_skyline:{
01760 if (gm.dsky==NULL){
01761 gm.dsky = new dskyline ();
01762 }
01763 dsky -> copy_dsky (*gm.dsky);
01764 break;
01765 }
01766 case compressed_rows:{
01767 if (gm.cr==NULL){
01768 gm.cr = new comprow ();
01769 }
01770 cr -> copy_cr (*gm.cr);
01771 break;
01772 }
01773 case symm_comp_rows:{
01774 if (gm.scr==NULL){
01775 gm.scr = new symcomprow ();
01776 }
01777 scr -> copy_scr (*gm.scr);
01778 break;
01779 }
01780 default:{
01781 print_err("unknown type of matrix storage is required",__FILE__,__LINE__,__func__);
01782 }
01783 }
01784 }
01785
01786 void gmatrix::a12block (gtopology *top,double *block,long nrdof,FILE *out)
01787 {
01788 if (pmat==0)
01789 prepmat2 (top,out);
01790
01791 switch (ts){
01792 case skyline_matrix:{
01793 sky->ldl_a12block (block,nrdof);
01794 break;
01795 }
01796 case spdirldl:{
01797 sdirect->GetA12block(block);
01798 break;
01799 }
01800 default:{
01801 print_err("wrong type of matrix storage is required",__FILE__,__LINE__,__func__);
01802 }
01803 }
01804
01805 }
01806
01807
01808
01809
01810
01811
01812
01813
01814 void gmatrix::printmat (FILE *out)
01815 {
01816 switch (ts){
01817 case dense_matrix:{
01818 dm->printmat (out);
01819 break;
01820 }
01821 case skyline_matrix:{
01822 sky->printmat (out);
01823 break;
01824 }
01825 case double_skyline:{
01826 dsky->printmat (out);
01827 break;
01828 }
01829 case compressed_rows:{
01830 cr->printmat (out);
01831 break;
01832 }
01833 case symm_comp_rows:{
01834 scr->printmat (out);
01835 break;
01836 }
01837 case element_matrices:{
01838 em->printmat (out);
01839 break;
01840 }
01841 case spdirect_stor_scr:{
01842 scr->printmat (out);
01843 break;
01844 }
01845 case spdirect_stor_cr:{
01846 cr->printmat (out);
01847 break;
01848 }
01849 default:{
01850 print_err("unknown type of matrix storage is required",__FILE__,__LINE__,__func__);
01851 }
01852 }
01853 }
01854
01855
01856
01857
01858
01859
01860
01861
01862 void gmatrix::printdiag (FILE *out)
01863 {
01864 switch (ts){
01865 case dense_matrix:{
01866 dm->printdiag (out);
01867 break;
01868 }
01869 case skyline_matrix:{
01870 sky->printdiag (out);
01871 break;
01872 }
01873 case double_skyline:{
01874 dsky->printdiag (out);
01875 break;
01876 }
01877 case compressed_rows:{
01878 cr->printdiag (out);
01879 break;
01880 }
01881 case symm_comp_rows:{
01882 scr->printdiag (out);
01883 break;
01884 }
01885 case spdirect_stor_cr:{
01886 cr->printdiag (out);
01887 break;
01888 }
01889 case spdirect_stor_scr:{
01890 scr->printdiag (out);
01891 break;
01892 }
01893 default:{
01894 print_err("unknown type of matrix storage is required",__FILE__,__LINE__,__func__);
01895 }
01896 }
01897 }
01898
01899 void gmatrix::changedecomp ()
01900 {
01901 switch (ts){
01902 case dense_matrix:{
01903 dm->changedecomp ();
01904 break;
01905 }
01906 case skyline_matrix:{
01907 sky->changedecomp ();
01908 break;
01909 }
01910 case double_skyline:{
01911 dsky->changedecomp ();
01912 break;
01913 }
01914 case compressed_rows:{
01915
01916 break;
01917 }
01918 case symm_comp_rows:{
01919
01920 break;
01921 }
01922 default:{
01923 print_err("unknown type of matrix storage is required",__FILE__,__LINE__,__func__);
01924 }
01925 }
01926
01927 }
01928
01929
01930
01931
01932
01933
01934
01935
01936
01937
01938
01939
01940
01941
01942
01943
01944
01945
01946
01947
01948
01949
01950
01951
01952
01953
01954
01955
01956
01957
01958
01959
01960
01961
01962
01963
01964
01965
01966
01967
01968
01969
01970
01971
01972
01973
01974
01975
01976
01977
01978
01979
01980
01981
01982
01983
01984
01985
01986
01987 long gmatrix::decomp ()
01988 {
01989 long d=0;
01990
01991 switch (ts){
01992 case dense_matrix:{
01993 d=dm->decomp ();
01994 break;
01995 }
01996 case skyline_matrix:{
01997 d=sky->decomp ();
01998 break;
01999 }
02000 case double_skyline:{
02001 d=dsky->decomp ();
02002 break;
02003 }
02004 default:{
02005 print_err("unknown type of matrix storage is required",__FILE__,__LINE__,__func__);
02006 }
02007 }
02008 return d;
02009 }
02010
02011
02012
02013
02014
02015
02016
02017
02018 double gmatrix::give_entry (long ri,long ci)
02019 {
02020 double e;
02021
02022 switch (ts){
02023 case dense_matrix:{
02024 e=dm->give_entry (ri,ci);
02025 break;
02026 }
02027 case skyline_matrix:{
02028 e=sky->give_entry (ri,ci);
02029 break;
02030 }
02031 case double_skyline:{
02032
02033 break;
02034 }
02035 case compressed_rows:{
02036 e=cr->give_entry (ri,ci);
02037 break;
02038 }
02039 case symm_comp_rows:{
02040 e=scr->give_entry (ri,ci);
02041 break;
02042 }
02043 default:{
02044 print_err("unknown type of matrix storage is required",__FILE__,__LINE__,__func__);
02045 }
02046 }
02047
02048 return e;
02049 }
02050
02051
02052
02053
02054
02055
02056
02057
02058
02059 void gmatrix::add_entry (double e,long ri,long ci)
02060 {
02061 switch (ts){
02062 case dense_matrix:{
02063 dm->add_entry (e,ri,ci);
02064 break;
02065 }
02066 case skyline_matrix:{
02067 sky->add_entry (e,ri,ci);
02068 break;
02069 }
02070 case double_skyline:{
02071
02072 break;
02073 }
02074 case compressed_rows:{
02075 cr->add_entry (e,ri,ci);
02076 break;
02077 }
02078 case symm_comp_rows:{
02079 scr->add_entry (e,ri,ci);
02080 break;
02081 }
02082 default:{
02083 print_err("unknown type of matrix storage is required",__FILE__,__LINE__,__func__);
02084 }
02085 }
02086
02087 }
02088
02089
02090
02091
02092
02093
02094 long gmatrix::give_negm ()
02095 {
02096 long negm;
02097
02098 switch (ts){
02099 case dense_matrix:{
02100 negm=dm->give_negm ();
02101 break;
02102 }
02103 case skyline_matrix:{
02104 negm=sky->give_negm ();
02105 break;
02106 }
02107 case double_skyline:{
02108 negm=dsky->give_negm ();
02109 break;
02110 }
02111 case compressed_rows:{
02112 negm=cr->give_negm ();
02113 break;
02114 }
02115 case symm_comp_rows:{
02116 negm=scr->give_negm ();
02117 break;
02118 }
02119 case spdirect_stor_cr:{
02120 negm=cr->give_negm ();
02121 break;
02122 }
02123 case spdirect_stor_scr:{
02124 negm=scr->give_negm ();
02125 break;
02126 }
02127 default:{
02128 print_err("unknown type of matrix storage is required",__FILE__,__LINE__,__func__);
02129 }
02130 }
02131
02132 return negm;
02133 }
02134
02135
02136
02137
02138
02139
02140 void gmatrix::diag_scale (double *d)
02141 {
02142 switch (ts){
02143 case dense_matrix:{
02144 dm->diag_scale (d);
02145 break;
02146 }
02147 case skyline_matrix:{
02148 sky->diag_scale (d);
02149 break;
02150 }
02151 default:{
02152 print_err("unknown type of matrix storage is required",__FILE__,__LINE__,__func__);
02153 }
02154 }
02155
02156 }
02157
02158
02159
02160
02161
02162
02163
02164
02165
02166
02167
02168 void gmatrix::diag_check (double thr,double *rhs)
02169 {
02170 switch (ts){
02171 case skyline_matrix:{
02172 sky->diag_check (thr);
02173 break;
02174 }
02175 case double_skyline:{
02176 dsky->diag_check (thr,rhs);
02177 break;
02178 }
02179 default:{
02180 print_err("unknown type of matrix storage is required",__FILE__,__LINE__,__func__);
02181 }
02182 }
02183
02184 }
02185
02186
02187
02188
02189
02190
02191
02192
02193
02194
02195
02196
02197 double gmatrix::power_method (double *v,long ni,double err)
02198 {
02199 long i;
02200 double ev,evo,nom,denom;
02201 double *p;
02202
02203 p = new double [n];
02204
02205 for (i=0;i<n;i++){
02206 v[i]=1.0;
02207 }
02208
02209 evo=0.0;
02210
02211 for (i=0;i<ni;i++){
02212 gmxv (v,p);
02213
02214
02215 nom=ss(v,p,n);
02216
02217 denom=ss(v,v,n);
02218
02219 ev=nom/denom;
02220
02221 fprintf (stdout,"\n power method, step %6ld, quotient %le",i,ev);
02222
02223 copyv(p,v,n);
02224
02225 cmulv(1.0/sqrt(denom),v,n);
02226
02227 if (fabs(ev-evo)/ev<err && i>0)
02228 break;
02229 evo=ev;
02230 }
02231
02232 delete [] p;
02233
02234 return ev;
02235 }
02236
02237
02238
02239
02240
02241
02242
02243
02244
02245
02246
02247 double gmatrix::inverse_iteration (double *v,long ni,double err)
02248 {
02249 long i;
02250 double ev,evo,nom,denom;
02251 double *p;
02252
02253 p = new double [n];
02254
02255 for (i=0;i<n;i++){
02256 v[i]=1.0;
02257 }
02258
02259 evo=0.0;
02260
02261 dm->lu (p,v,zero,2);
02262
02263 for (i=0;i<ni;i++){
02264 dm->lu (p,v,zero,3);
02265
02266
02267 nom=ss(v,p,n);
02268
02269 denom=ss(v,v,n);
02270
02271 ev=nom/denom;
02272
02273 fprintf (stdout,"\n inverse iteration, step %6ld, quotient %le",i,ev);
02274
02275 copyv(p,v,n);
02276
02277 cmulv(1.0/sqrt(denom),v,n);
02278
02279 if (fabs(ev-evo)/ev<err && i>0){
02280 break;
02281 }
02282
02283 evo=ev;
02284 }
02285
02286 delete [] p;
02287
02288 return ev;
02289
02290 }
02291
02292
02293
02294
02295
02296
02297
02298
02299
02300
02301
02302 double gmatrix::estim_spect_radius ()
02303 {
02304 double sr=0.0;
02305
02306 switch (ts){
02307 case dense_matrix:{
02308 sr=dm->estim_spect_radius ();
02309 break;
02310 }
02311 case compressed_rows:{
02312 sr=cr->estim_spect_radius ();
02313 break;
02314 }
02315 default:{
02316 print_err("unknown type of matrix storage is required",__FILE__,__LINE__,__func__);
02317 }
02318 }
02319
02320 return sr;
02321 }
02322