00001 #include "densemat.h"
00002 #include <math.h>
00003 #include <string.h>
00004 #include <stdlib.h>
00005
00006 #include "precond.h"
00007
00008 densemat::densemat (void)
00009 {
00010
00011 n=0;
00012
00013 a=NULL;
00014
00015 decompid=0;
00016 }
00017
00018 densemat::~densemat (void)
00019 {
00020 delete [] a;
00021 }
00022
00023
00024
00025
00026
00027
00028
00029
00030 void densemat::alloc (long m)
00031 {
00032 n=m;
00033 negm=m*m;
00034 a = new double [negm];
00035 memset (a,0,negm*sizeof(double));
00036 }
00037
00038 void densemat::dealloc (void)
00039 {
00040 delete [] a;
00041 a = NULL;
00042 }
00043
00044 void densemat::copy (densemat *dm)
00045 {
00046 long i,j;
00047
00048 if (a != NULL)
00049 delete [] a;
00050
00051 n=dm->n;
00052 negm=n*n;
00053 a = new double [negm];
00054
00055 for (i=0;i<n;i++){
00056 for (j=0;j<n;j++){
00057 a[i*n+j]=dm->a[i*n+j];
00058 }
00059 }
00060 }
00061
00062
00063
00064
00065
00066
00067
00068
00069 void densemat::copy_dm (densemat &dm)
00070 {
00071 long i,j;
00072
00073
00074 dm.n=n;
00075
00076 dm.negm=n*n;
00077
00078 if (dm.a != NULL)
00079 delete [] dm.a;
00080 dm.a = new double [negm];
00081
00082 for (i=0;i<n;i++){
00083 for (j=0;j<n;j++){
00084 dm.a[i*n+j]=a[i*n+j];
00085 }
00086 }
00087 }
00088
00089
00090
00091
00092
00093
00094
00095 double* densemat::status ()
00096 {
00097 return a;
00098 }
00099
00100
00101
00102
00103
00104
00105 long densemat::decomp ()
00106 {
00107 return decompid;
00108 }
00109
00110
00111
00112
00113
00114
00115 void densemat::changedecomp ()
00116 {
00117 if (decompid==0) decompid=1;
00118 else decompid=0;
00119 }
00120
00121 void densemat::setfact ()
00122 {
00123 decompid=1;
00124 }
00125 void densemat::setnotfact ()
00126 {
00127 decompid=0;
00128 }
00129
00130
00131
00132
00133
00134
00135 void densemat::nullmat ()
00136 {
00137 long i,j;
00138 for (i=0;i<n;i++){
00139 for (j=0;j<n;j++){
00140 a[i*n+j]=0.0;
00141 }
00142 }
00143 }
00144
00145
00146
00147
00148
00149
00150
00151
00152
00153 void densemat::localize (matrix &b,long *cn)
00154 {
00155 long i,j,ii,jj;
00156
00157 for (i=0;i<b.m;i++){
00158 ii=cn[i]-1;
00159 if (ii<0) continue;
00160 for (j=0;j<b.m;j++){
00161 jj=cn[j]-1;
00162 if (jj<0) continue;
00163 a[ii*n+jj]+=b[i][j];
00164 }
00165 }
00166 }
00167
00168
00169
00170
00171
00172
00173
00174
00175
00176
00177 void densemat::localized (double *b,long *cn,long m)
00178 {
00179 long i,j,ii,jj;
00180
00181 for (i=0;i<m;i++){
00182 ii=cn[i]-1;
00183 if (ii<0) continue;
00184 for (j=0;j<m;j++){
00185 jj=cn[j]-1;
00186 if (jj<0) continue;
00187 a[ii*n+jj]+=b[i*m+j];
00188 }
00189 }
00190 }
00191
00192
00193
00194
00195
00196
00197
00198
00199
00200
00201 void densemat::glocalize (matrix &b,long *rcn,long *ccn)
00202 {
00203 long i,j,ii,jj;
00204
00205 for (i=0;i<b.m;i++){
00206 ii=rcn[i]-1;
00207 if (ii<0) continue;
00208 for (j=0;j<b.n;j++){
00209 jj=ccn[j]-1;
00210 if (jj<0) continue;
00211 a[ii*n+jj]+=b[i][j];
00212 }
00213 }
00214 }
00215
00216
00217
00218
00219
00220
00221
00222
00223
00224
00225
00226 void densemat::mult_localize (long nm,long *ncn1,long *ncn2,long *mcn)
00227 {
00228 long i,cn,m;
00229
00230 for (i=0;i<nm;i++){
00231 cn=mcn[i]-1;
00232 if (cn>-1){
00233 m=ncn1[i]-1;
00234 if (m>-1){
00235 a[m*n+cn]=-1.0;
00236 a[cn*n+m]=-1.0;
00237 }
00238 m=ncn2[i]-1;
00239 if (m>-1){
00240 a[m*n+cn]=1.0;
00241 a[cn*n+m]=1.0;
00242 }
00243 }
00244 }
00245 }
00246
00247
00248
00249
00250
00251
00252
00253
00254
00255 void densemat::initiate (long ndof,long mespr)
00256 {
00257 if (status ()==NULL){
00258 alloc (ndof);
00259 }
00260 else{
00261 nullmat ();
00262 }
00263
00264 if (mespr==1) fprintf (stdout,"\n number of matrix entries %ld",negm);
00265 }
00266
00267
00268
00269
00270
00271
00272
00273
00274
00275
00276 void densemat::mxv_dm (double *b,double *c)
00277 {
00278 long i,j,aca;
00279 double s;
00280
00281 aca=0;
00282 for (i=0;i<n;i++){
00283 s=0.0;
00284 for (j=0;j<n;j++){
00285 s+=a[aca]*b[j];
00286 aca++;
00287 }
00288 c[i]=s;
00289 }
00290 }
00291
00292
00293
00294
00295
00296
00297
00298
00299
00300 void densemat::addmat_dm (double c,densemat &dm)
00301 {
00302 long i,j;
00303 for (i=0;i<n;i++){
00304 for (j=0;j<n;j++){
00305 a[i*n+j]+=c*dm.a[i*n+j];
00306 }
00307 }
00308 }
00309
00310
00311
00312
00313
00314
00315
00316
00317 void densemat::scalmat_dm (double c)
00318 {
00319 long i,j;
00320 for (i=0;i<n;i++){
00321 for (j=0;j<n;j++){
00322 a[i*n+j]*=c;
00323 }
00324 }
00325 }
00326
00327
00328
00329
00330
00331
00332
00333
00334
00335
00336 void densemat::maxmb(long nm, double *ma, double *mb, double *mc)
00337 {
00338 long i,j,k;
00339 double s;
00340
00341 for(i = 0; i < nm; i++){
00342 for(j = 0; j < nm; j++){
00343 s=0.0;
00344 for(k = 0; k< nm; k++){
00345 s+=ma[i*nm+k]*mb[k*nm+j];
00346 }
00347 mc[i*n+j]=s;
00348 }
00349 }
00350
00351
00352 }
00353
00354
00355
00356
00357
00358
00359
00360
00361
00362
00363
00364
00365
00366
00367
00368
00369
00370
00371
00372 void densemat::gemp (double *x,double *y,long m,double zero,long pivot)
00373 {
00374 long i,j,k,ii,jj,kk,ir,ic,li,ui,*order;
00375 double f,s;
00376
00377 order = new long [n];
00378
00379
00380
00381 for (i=0;i<n;i++){
00382 order[i]=i;
00383 }
00384
00385 for (k=0;k<n-1;k++){
00386 ir=k; ic=k;
00387
00388 if (pivot==1){
00389
00390 if (fabs(a[k*n+k])<zero){
00391 s=0.0; ii=k*n+k;
00392 for (i=k;i<n;i++){
00393 for (j=k;j<n;j++){
00394 f=fabs(a[ii]);
00395 if (f>s){
00396 s=f; ir=i; ic=j;
00397 }
00398 ii++;
00399 }
00400 ii+=k;
00401 }
00402 if (s<zero){
00403 print_err("singular matrix is detected", __FILE__, __LINE__, __func__);
00404 break;
00405 }
00406 }
00407 }
00408
00409 if (pivot==2){
00410
00411 s=0.0; ii=k*n+k;
00412 for (i=k;i<n;i++){
00413 for (j=k;j<n;j++){
00414 f=fabs(a[ii]);
00415 if (f>s){
00416 s=f; ir=i; ic=j;
00417 }
00418 ii++;
00419 }
00420 ii+=k;
00421 }
00422 if (s<zero){
00423 print_err("singular matrix is detected", __FILE__, __LINE__, __func__);
00424 break;
00425 }
00426 }
00427
00428
00429 if (ir!=k){
00430 li=k*n+k; ui=(k+1)*n; j=ir*n+k;
00431 for (i=li;i<ui;i++){
00432 s=a[i];
00433 a[i]=a[j];
00434 a[j]=s;
00435 j++;
00436 }
00437
00438 li=k*m; ui=li+m; j=ir*m;
00439 for (i=li;i<ui;i++){
00440 s=y[i];
00441 y[i]=y[j];
00442 y[j]=s;
00443 j++;
00444 }
00445 }
00446
00447
00448 if (ic!=k){
00449 i=order[k];
00450 order[k]=order[ic];
00451 order[ic]=i;
00452
00453 j=k; ii=ic;
00454 for (i=0;i<n;i++){
00455 s=a[j];
00456 a[j]=a[ii];
00457 a[ii]=s;
00458 j+=n; ii+=n;
00459 }
00460 }
00461
00462
00463
00464
00465 for (i=k+1;i<n;i++){
00466 ii=i*n+k; kk=k*n+k;
00467 s=a[ii]/a[kk];
00468
00469
00470 for (j=k;j<n;j++){
00471 a[ii]-=s*a[kk];
00472 ii++; kk++;
00473 }
00474
00475
00476 ii=i*m; kk=k*m;
00477 for (j=0;j<m;j++){
00478 y[ii]-=s*y[kk];
00479 ii++; kk++;
00480 }
00481 }
00482
00483
00484 }
00485
00486 if (fabs(a[(n-1)*n+n-1])<zero){
00487 print_err("singular matrix is detected", __FILE__, __LINE__, __func__);
00488 }
00489
00490
00491
00492
00493 for (k=n-1;k>-1;k--){
00494 f=a[k*n+k]; kk=k*m;
00495 for (i=0;i<m;i++){
00496 s=0.0; ii=k*n+n-1; jj=(n-1)*m+i;
00497 for (j=n-1;j>k;j--){
00498 s+=a[ii]*x[jj];
00499 ii--; jj-=m;
00500 }
00501 x[kk]=(y[kk]-s)/f;
00502 kk++;
00503 }
00504 }
00505
00506
00507 for (i=0;i<n;i++){
00508 if (order[i]!=i){
00509 for (j=i;j<n;j++){
00510 if (order[j]==i){
00511 jj=j; break;
00512 }
00513 }
00514
00515 ii=order[i];
00516 order[i]=order[jj];
00517 order[jj]=ii;
00518
00519 ii=i*m; kk=jj*m;
00520 for (j=0;j<m;j++){
00521 s=x[ii];
00522 x[ii]=x[kk];
00523 x[kk]=s;
00524 ii++; kk++;
00525 }
00526 }
00527 }
00528
00529 delete [] order;
00530 }
00531
00532
00533
00534
00535
00536
00537
00538
00539
00540
00541
00542
00543
00544
00545
00546
00547
00548
00549
00550
00551
00552
00553 void densemat::gempkon (double *b,double *c,double *x,double *y,long m,double zero,long tc)
00554 {
00555 long i,j,k,ii,kk;
00556 double s;
00557
00558 if (tc==1){
00559 for (k=0;k<n-m;k++){
00560
00561 if (fabs(a[k*n+k])<zero) fprintf (stderr,"\n");
00562
00563 for (i=k+1;i<n;i++){
00564 ii=i*n+k; kk=k*n+k;
00565 s=a[ii]/a[kk];
00566
00567
00568 for (j=k;j<n;j++){
00569 a[ii]-=s*a[kk];
00570 ii++; kk++;
00571 }
00572
00573
00574 y[i]-=s*y[k];
00575
00576 }
00577
00578 }
00579
00580
00581 ii=n-m;
00582 for (i=0;i<m;i++){
00583 kk=n-m;
00584 for (j=0;j<m;j++){
00585 b[i*m+j]=a[ii*n+kk];
00586 kk++;
00587 }
00588 ii++;
00589 }
00590
00591
00592 ii=n-m;
00593 for (i=0;i<m;i++){
00594 c[i]=y[ii];
00595 ii++;
00596 }
00597
00598
00599 }
00600
00601 if (tc==2){
00602
00603 j=0;
00604 for (k=n-m;k<n;k++){
00605 x[k]=c[j]; j++;
00606 }
00607
00608
00609 for (k=n-m-1;k>-1;k--){
00610 s=0.0; ii=k*n+k+1;
00611 for (j=k+1;j<n;j++){
00612 s+=a[ii]*x[j];
00613 ii++;
00614 }
00615 x[k]=(y[k]-s)/a[k*n+k];
00616 }
00617 }
00618
00619 }
00620
00621
00622
00623
00624
00625
00626
00627
00628
00629
00630
00631
00632
00633
00634
00635 void densemat::lu (double *x,double *y,double zero,long tc)
00636 {
00637 long i,j,k;
00638 double s;
00639
00640 if (tc==1 || tc==2){
00641 if (decompid==0){
00642 for (i=0;i<n;i++){
00643 for (j=0;j<i;j++){
00644 s=0.0;
00645 for (k=0;k<j;k++){
00646 s+=a[i*n+k]*a[k*n+j];
00647 }
00648 a[i*n+j]-=s;
00649 }
00650
00651 s=0.0;
00652 for (k=0;k<i;k++){
00653 s+=a[i*n+k]*a[k*n+i];
00654 }
00655 a[i*n+i]-=s;
00656 if (fabs(a[i*n+i])<zero){
00657 print_err("zero diagonal entry in lu factorization", __FILE__, __LINE__, __func__);
00658 }
00659
00660 for (j=i+1;j<n;j++){
00661 s=0.0;
00662 for (k=0;k<i;k++){
00663 s+=a[i*n+k]*a[k*n+j];
00664 }
00665 a[i*n+j]=(a[i*n+j]-s)/a[i*n+i];
00666 }
00667 }
00668 decompid=1;
00669 }
00670 }
00671
00672 if (tc==1 || tc==3){
00673 if (decompid==0){
00674 print_err("matrix is not factorized, back substitution cannot be performed", __FILE__, __LINE__, __func__);
00675 }
00676 else{
00677 for (i=0;i<n;i++){
00678 s=0.0;
00679 for (j=0;j<i;j++){
00680 s+=a[i*n+j]*y[j];
00681 }
00682 y[i]=(y[i]-s)/a[i*n+i];
00683 }
00684 for (i=n-1;i>-1;i--){
00685 s=0.0;
00686 for (j=i+1;j<n;j++){
00687 s+=a[i*n+j]*x[j];
00688 }
00689 x[i]=y[i]-s;
00690 }
00691 }
00692 }
00693 }
00694
00695
00696
00697
00698
00699
00700
00701
00702
00703
00704
00705
00706
00707
00708
00709
00710
00711
00712
00713
00714 void densemat::ll (double *x,double *y,double zero,long tc)
00715 {
00716 long i,j,k;
00717 double s;
00718 double *z;
00719
00720 z = new double [n];
00721
00722 if (tc==1 || tc==2){
00723 if (decompid==0){
00724 for (i=0;i<n;i++){
00725 for (j=0;j<i;j++){
00726 s=0.0;
00727 for (k=0;k<j;k++){
00728 s+=a[i*n+k]*a[j*n+k];
00729 }
00730 a[i*n+j]-=s;
00731 a[i*n+j]/=a[j*n+j];
00732 }
00733
00734 s=0.0;
00735 for (k=0;k<i;k++){
00736 s+=a[i*n+k]*a[i*n+k];
00737 }
00738 a[i*n+i]-=s;
00739 if (a[i*n+i]<0.0){
00740 print_err("negative diagonal entry in Cholesky factorization", __FILE__, __LINE__, __func__);
00741 }
00742 a[i*n+i]=sqrt(a[i*n+i]);
00743
00744 }
00745
00746 decompid=1;
00747 }
00748 }
00749
00750 if (tc==1 || tc==3){
00751 if (decompid==0){
00752 print_err("matrix is not factorized, back substitution cannot be performed", __FILE__, __LINE__, __func__);
00753 }
00754 else{
00755 for (i=0;i<n;i++){
00756 s=0.0;
00757 for (k=0;k<i;k++){
00758 s+=a[i*n+k]*z[k];
00759 }
00760 z[i]=(y[i]-s)/a[i*n+i];
00761 }
00762
00763 for (i=n-1;i>-1;i--){
00764 s=0.0;
00765 for (k=i+1;k<n;k++){
00766 s+=a[k*n+i]*x[k];
00767 }
00768 x[i]=(z[i]-s)/a[i*n+i];
00769 }
00770 }
00771 }
00772
00773 delete [] z;
00774 }
00775
00776
00777
00778
00779
00780
00781
00782
00783
00784
00785
00786
00787
00788
00789
00790
00791
00792
00793
00794
00795
00796 void densemat::ill (double *x,double *y,double zero,double limit,long tc)
00797 {
00798 long i,j,k;
00799 double s;
00800 double *z;
00801
00802 z = new double [n];
00803
00804 if (tc==1 || tc==2){
00805 if (decompid==0){
00806 for (i=0;i<n;i++){
00807 for (j=0;j<i;j++){
00808 if (fabs(a[i*n+j])>limit){
00809 s=0.0;
00810 for (k=0;k<j;k++){
00811 s+=a[i*n+k]*a[j*n+k];
00812 }
00813 a[i*n+j]-=s;
00814 a[i*n+j]/=a[j*n+j];
00815 }
00816 }
00817
00818 s=0.0;
00819 for (k=0;k<i;k++){
00820 s+=a[i*n+k]*a[i*n+k];
00821 }
00822 a[i*n+i]-=s;
00823 if (a[i*n+i]<0.0){
00824 print_err("negative diagonal entry in Cholesky decomposition", __FILE__, __LINE__, __func__);
00825 }
00826 a[i*n+i]=sqrt(a[i*n+i]);
00827
00828 }
00829 decompid=1;
00830 }
00831 }
00832
00833 if (tc==1 || tc==3){
00834 if (decompid==0){
00835 print_err("matrix is not factorized, back substitution cannot be performed", __FILE__, __LINE__, __func__);
00836 }
00837 else{
00838 for (i=0;i<n;i++){
00839 s=0.0;
00840 for (k=0;k<i;k++){
00841 s+=a[i*n+k]*z[k];
00842 }
00843 z[i]=(y[i]-s)/a[i*n+i];
00844 }
00845
00846 for (i=n-1;i>-1;i--){
00847 s=0.0;
00848 for (k=i+1;k<n;k++){
00849 s+=a[k*n+i]*x[k];
00850 }
00851 x[i]=(z[i]-s)/a[i*n+i];
00852 }
00853 }
00854 }
00855
00856 delete [] z;
00857 }
00858
00859
00860
00861
00862
00863
00864
00865
00866
00867
00868
00869
00870
00871
00872 void densemat::ker (double *r,long &dim,long *se,long ense,double limit)
00873 {
00874 long i,j,k,stop;
00875 double s,*b;
00876
00877
00878 b = new double [ense*n];
00879
00880
00881 dim=0;
00882 for (i=0;i<n;i++){
00883 if (fabs(a[i*n+i])<limit){
00884 se[dim]=i;
00885
00886 for (j=0;j<n;j++){
00887 a[i*n+j]=0.0;
00888 }
00889
00890 a[i*n+i]=1.0;
00891 k=dim*n;
00892 for (j=0;j<n;j++){
00893 b[k]=a[j*n+i]; k++;
00894 a[j*n+i]=0.0;
00895 }
00896 a[i*n+i]=1.0;
00897
00898 dim++;
00899 if (dim>ense){
00900 print_err("wrong number of singular equations, there are more singular equations than has been estimated", __FILE__, __LINE__, __func__);
00901 abort ();
00902 }
00903 }
00904 else{
00905 for (j=i+1;j<n;j++){
00906 s=a[j*n+i]/a[i*n+i];
00907 for (k=i;k<n;k++){
00908 a[j*n+k]-=s*a[i*n+k];
00909 }
00910 }
00911 }
00912 }
00913
00914
00915
00916 k=0;
00917 for (i=0;i<dim;i++){
00918 for (j=0;j<n;j++){
00919 a[j*n+se[i]]=b[k]; k++;
00920 }
00921 }
00922
00923 nullv (r,dim*n);
00924 for (i=0;i<dim;i++){
00925 r[i*n+se[i]]=1.0;
00926 for (j=n-1;j>-1;j--){
00927 stop=0;
00928 for (k=0;k<dim;k++){
00929 if (j==se[k]){
00930 stop=1; break;
00931 }
00932 }
00933 if (stop==1) continue;
00934
00935 s=0.0;
00936 for (k=j+1;k<n;k++){
00937 s-=a[j*n+k]*r[i*n+k];
00938 }
00939
00940 r[i*n+j]=s/a[j*n+j];
00941 }
00942 }
00943 delete [] b;
00944
00945 }
00946
00947
00948
00949
00950
00951
00952
00953
00954
00955
00956
00957
00958
00959
00960
00961
00962
00963
00964
00965
00966
00967
00968
00969 void densemat::cg (double *x,double *y,long ni,double err,long &ani,double &ares,double zero,long iv)
00970 {
00971 long i,j;
00972 double nom,denom,nory,alpha,beta;
00973 double *d,*r,*p;
00974
00975 d = new double [n];
00976 r = new double [n];
00977 p = new double [n];
00978
00979
00980 if (iv==0){
00981 for (i=0;i<n;i++)
00982 x[i]=0.0;
00983 }
00984 mxv_dm (x,p);
00985
00986 nory=0.0; nom=0.0;
00987 for (i=0;i<n;i++){
00988 nory+=y[i]*y[i];
00989 r[i]=p[i]-y[i];
00990 nom+=r[i]*r[i];
00991 d[i]=-1.0*r[i];
00992 }
00993
00994 if (nory<zero){
00995 print_err("norm of right hand side in conjugate gradient method is smaller than %e", __FILE__, __LINE__, __func__,zero);
00996 ares=nory; ani=0;
00997 for (i=0;i<n;i++)
00998 x[i]=0.0;
00999 return;
01000 }
01001
01002
01003
01004 for (i=0;i<ni;i++){
01005
01006
01007
01008 mxv_dm (d,p);
01009
01010
01011 denom = ss (d,p,n);
01012 if (fabs(denom)<zero){
01013 fprintf (stdout,"\n there is zero denominator in alpha computation in conjugate gradient method (cr.cpp)\n");
01014 break;
01015 }
01016
01017 alpha = nom/denom;
01018
01019
01020 for (j=0;j<n;j++){
01021 x[j]+=alpha*d[j];
01022 r[j]+=alpha*p[j];
01023 }
01024
01025 denom=nom;
01026
01027 nom = ss (r,r,n);
01028
01029 fprintf (stdout,"\n iteration %ld norres/norrhs %e",i,nom/nory);
01030
01031 if (nom/nory<err) break;
01032
01033
01034
01035 beta = nom/denom;
01036
01037
01038 for (j=0;j<n;j++){
01039 d[j]=beta*d[j]-r[j];
01040 }
01041 }
01042
01043 ani=i; ares=nom;
01044
01045 delete [] p; delete [] r; delete [] d;
01046
01047 }
01048
01049
01050
01051
01052
01053
01054
01055
01056
01057
01058
01059
01060
01061
01062
01063
01064
01065
01066
01067
01068
01069
01070
01071
01072 void densemat::cg_prec (densemat &dm,double *x,double *y,long ni,double err,long &ani,double &ares,
01073 double zero,long iv,long tprec,double par)
01074 {
01075 long i,j;
01076 double nom,denom,nory,alpha,beta;
01077 double *d,*r,*p,*h;
01078
01079 d = new double [n];
01080 r = new double [n];
01081 p = new double [n];
01082 h = new double [n];
01083
01084
01085 if (iv==0){
01086 for (i=0;i<n;i++)
01087 x[i]=0.0;
01088 }
01089 mxv_dm (x,p);
01090 nory=0.0;
01091 for (i=0;i<n;i++){
01092 nory+=y[i]*y[i];
01093 r[i]=p[i]-y[i];
01094 }
01095
01096 if (nory<zero){
01097 print_err("right hand side vector has zero norm", __FILE__, __LINE__, __func__);
01098 ares=nory; ani=0;
01099 }
01100
01101 switch (tprec){
01102
01103
01104 case 10:{
01105
01106 dm.ill (h,r,zero,zero,3);
01107 break;
01108 }
01109
01110
01111
01112
01113 default:{
01114 print_err("wrong preconditioner type is required", __FILE__, __LINE__, __func__);
01115 }
01116 }
01117
01118 nom = ss (r,h,n);
01119
01120 for (i=0;i<n;i++){
01121 d[i]=0.0-h[i];
01122 }
01123
01124
01125 for (i=0;i<ni;i++){
01126
01127
01128 mxv_dm (d,p);
01129
01130 denom = ss (d,p,n);
01131 if (fabs(denom)<zero){
01132 print_err("denominator in alpha expression is equal to zero", __FILE__, __LINE__, __func__);
01133 ares=nom; ani=i;
01134 return;
01135 }
01136
01137
01138 alpha = nom/denom;
01139
01140
01141 for (j=0;j<n;j++){
01142 x[j]+=alpha*d[j];
01143 r[j]+=alpha*p[j];
01144 }
01145
01146
01147
01148 switch (tprec){
01149
01150
01151 case 10:{
01152
01153 dm.ill (h,r,zero,zero,3);
01154 break;
01155 }
01156
01157
01158
01159
01160 default:{
01161 print_err("wrong preconditioner type is required", __FILE__, __LINE__, __func__);
01162 }
01163 }
01164
01165 denom=nom;
01166 nom = ss (r,h,n);
01167 beta = nom/denom;
01168
01169 fprintf (stdout,"\n iteration %ld norres/norrhs %e",i,nom/nory);
01170
01171 if (fabs(nom/nory)<err){
01172 ani=i; ares=nom/nory;
01173 break;
01174 }
01175
01176
01177 for (j=0;j<n;j++){
01178 d[j]=beta*d[j]-h[j];
01179 }
01180 }
01181
01182 delete [] h;
01183 delete [] p;
01184 delete [] r;
01185 delete [] d;
01186 }
01187
01188
01189
01190
01191
01192
01193
01194
01195 void densemat::cg_prec_new (precond &pr,double *x,double *y,long ni,double err,long &ani,double &ares,
01196 double zero,long iv)
01197 {
01198 long i,j;
01199 double nom,denom,nory,alpha,beta;
01200 double *d,*r,*p,*h;
01201
01202 d = new double [n];
01203 r = new double [n];
01204 p = new double [n];
01205 h = new double [n];
01206
01207
01208 if (iv==0){
01209 for (i=0;i<n;i++)
01210 x[i]=0.0;
01211 }
01212
01213
01214 mxv_dm (x,p);
01215
01216
01217 nory=0.0;
01218 for (i=0;i<n;i++){
01219 nory+=y[i]*y[i];
01220 r[i]=y[i]-p[i];
01221 }
01222
01223 if (nory<zero){
01224 print_err("right hand side vector has zero norm", __FILE__, __LINE__, __func__);
01225 ares=nory; ani=0;
01226 }
01227
01228
01229 pr.preconditioning (r,h,y);
01230
01231 nom = ss (r,h,n);
01232
01233
01234 for (i=0;i<n;i++){
01235 d[i]=h[i];
01236 }
01237
01238
01239
01240 for (i=0;i<ni;i++){
01241
01242
01243 mxv_dm (d,p);
01244
01245 denom = ss (d,p,n);
01246 if (fabs(denom)<zero){
01247 print_err("denominator in alpha expression is equal to zero", __FILE__, __LINE__, __func__);
01248 ares=nom; ani=i;
01249 return;
01250 }
01251
01252
01253 alpha = nom/denom;
01254
01255
01256 for (j=0;j<n;j++){
01257 x[j]+=alpha*d[j];
01258 r[j]-=alpha*p[j];
01259 }
01260
01261
01262 pr.preconditioning (r,h,y);
01263
01264
01265 denom=nom;
01266 nom = ss (r,h,n);
01267 beta = nom/denom;
01268
01269 fprintf (stdout,"\n iteration %ld norres/norrhs %e",i,nom/nory);
01270
01271 if (fabs(nom/nory)<err){
01272 ani=i; ares=nom/nory;
01273 break;
01274 }
01275
01276
01277 for (j=0;j<n;j++){
01278 d[j]=beta*d[j]+h[j];
01279 }
01280 }
01281
01282 delete [] h;
01283 delete [] p;
01284 delete [] r;
01285 delete [] d;
01286 }
01287
01288
01289
01290
01291
01292
01293
01294
01295
01296
01297
01298
01299
01300
01301
01302
01303
01304
01305
01306 void densemat::printmat (FILE *out)
01307 {
01308 long i,j;
01309
01310
01311 fprintf (out,"\n\n");
01312 for (i=0;i<n;i++){
01313
01314 fprintf (out,"\n");
01315 for (j=0;j<n;j++){
01316
01317 fprintf (out," %20.15le",a[i*n+j]);
01318 }
01319 fprintf (out,";");
01320 }
01321
01322
01323
01324
01325
01326
01327
01328
01329
01330
01331
01332
01333
01334
01335
01336
01337
01338
01339 }
01340
01341
01342
01343
01344
01345
01346
01347
01348
01349 void densemat::printdiag (FILE *out)
01350 {
01351 long i;
01352
01353 fprintf (out,"\n\n");
01354 for (i=0;i<n;i++){
01355 fprintf (out,"%5ld %.10e\n",i,a[i*n+i]);
01356 }
01357 }
01358
01359
01360
01361
01362
01363
01364
01365
01366 double densemat::give_entry (long ri,long ci)
01367 {
01368 double e;
01369
01370 return e = a[ri*n+ci];
01371 }
01372
01373
01374
01375
01376
01377
01378
01379
01380
01381 void densemat::add_entry (double e,long ri,long ci)
01382 {
01383 a[ri*n+ci]+=e;
01384 }
01385
01386
01387
01388
01389
01390
01391 long densemat::give_negm ()
01392 {
01393 return negm;
01394 }
01395
01396
01397
01398
01399
01400
01401
01402
01403
01404 void densemat::diag_scale (double *d)
01405 {
01406 long i,j;
01407
01408 for (i=0;i<n;i++){
01409 d[i]=1.0/sqrt(a[i*n+i]);
01410 }
01411
01412 for (i=0;i<n;i++){
01413 for (j=0;j<n;j++){
01414 a[i*n+j]*=d[i]*d[j];
01415 }
01416 }
01417 }
01418
01419
01420
01421
01422
01423
01424
01425
01426
01427
01428 double densemat::estim_spect_radius ()
01429 {
01430 long i,j;
01431 double s,sr;
01432
01433
01434 sr=0.0;
01435 for (i=0;i<n;i++){
01436 s=0.0;
01437 for (j=0;j<n;j++){
01438 s+=fabs(a[i*n+j]);
01439 }
01440 if (s>sr)
01441 sr=s;
01442 }
01443
01444 return sr;
01445 }
01446
01447
01448
01449
01450
01451
01452
01453
01454
01455
01456 void densemat::assemble_dense_from_scr(long *scradr,long *scrci,double *scra,long neq)
01457 {
01458 long i,j;
01459 long address1,address2,index;
01460
01461
01462 n = neq;
01463
01464 a = new double [n*n];
01465
01466 for(i = 0; i < n*n; i++){
01467 a[i] = 0.0;
01468 }
01469
01470 for(i = 0; i < n; i++){
01471 address1 = scradr[i];
01472 address2 = scradr[i+1];
01473 for(j = address1; j < address2; j++){
01474 index = scrci[j];
01475 a[i*n+index]=scra[j];
01476 a[index*n+i]=scra[j];
01477 }
01478 }
01479
01480 }
01481