00001 #include "cr.h"
00002 #include <limits.h>
00003 #include <math.h>
00004 #include <string.h>
00005 #include <stdlib.h>
00006
00007 #include "precond.h"
00008
00009 comprow::comprow (void)
00010 {
00011
00012 n=0;
00013
00014 negm=0;
00015
00016 limit = 0.0;
00017
00018 decompid=0;
00019
00020 adr=NULL;
00021 ci=NULL;
00022 a=NULL;
00023 adra=NULL;
00024 aux=NULL;
00025 }
00026
00027
00028
00029 comprow::~comprow ( void )
00030 {
00031 delete [] adr;
00032 delete [] ci;
00033 delete [] a;
00034
00035 delete [] adra;
00036
00037 if (aux != NULL)
00038 delete [] aux;
00039 }
00040
00041
00042
00043
00044
00045
00046 long comprow::decomp ()
00047 {
00048 return decompid;
00049 }
00050
00051
00052
00053
00054
00055
00056 void comprow::changedecomp ()
00057 {
00058 if (decompid==0) decompid=1;
00059 else decompid=0;
00060 }
00061
00062
00063
00064
00065
00066
00067
00068
00069
00070 void comprow::allocadr (long m)
00071 {
00072 n=m;
00073 adr = new long [n+1];
00074 adra = new long [n+1];
00075
00076 memset (adr,0,(n+1)*sizeof(long));
00077 memset (adra,0,(n+1)*sizeof(long));
00078 }
00079
00080
00081
00082
00083
00084
00085 double* comprow::status ()
00086 {
00087 return a;
00088 }
00089
00090
00091
00092
00093
00094
00095 void comprow::nullmat ()
00096 {
00097 long i;
00098
00099 for (i=0;i<negm;i++){
00100 a[i]=0.0;
00101 }
00102 }
00103
00104
00105
00106
00107
00108
00109
00110
00111
00112 void comprow::numcontr_elem (long *cn,long ndofe)
00113 {
00114 long i,j,ii,jj;
00115
00116 for (i=0;i<ndofe;i++){
00117 ii=cn[i]-1;
00118 if (ii<0) continue;
00119 for (j=0;j<ndofe;j++){
00120 jj=cn[j]-1;
00121 if (jj<0) continue;
00122 adr[ii]++;
00123 }
00124 }
00125 }
00126
00127
00128
00129
00130
00131
00132
00133
00134 void comprow::numcontr (gtopology *top)
00135 {
00136 long i,ne,ndofe,*cn;
00137 ne=top->ne;
00138 for (i=0;i<ne;i++){
00139 ndofe=top->give_ndofe (i);
00140 cn = new long [ndofe];
00141 top->give_code_numbers (i,cn);
00142 numcontr_elem (cn,ndofe);
00143 delete [] cn;
00144 }
00145 }
00146
00147
00148
00149
00150
00151
00152
00153
00154
00155
00156
00157
00158 void comprow::fillarray_elem (long *cn,long ndofe)
00159 {
00160 long i,j,ii,jj;
00161
00162 for (i=0;i<ndofe;i++){
00163 ii=cn[i]-1;
00164 if (ii<0) continue;
00165 for (j=0;j<ndofe;j++){
00166 jj=cn[j]-1;
00167 if (jj<0) continue;
00168 aux[adra[ii]]=jj;
00169 adra[ii]++;
00170 }
00171 }
00172 }
00173
00174
00175
00176
00177
00178
00179
00180
00181
00182
00183
00184 void comprow::fillarray (gtopology *top)
00185 {
00186 long i,ne,ndofe,*cn;
00187 ne=top->ne;
00188 for (i=0;i<ne;i++){
00189 ndofe=top->give_ndofe (i);
00190 cn = new long [ndofe];
00191 top->give_code_numbers (i,cn);
00192 fillarray_elem (cn,ndofe);
00193 delete [] cn;
00194 }
00195 }
00196
00197
00198
00199
00200
00201
00202 void comprow::addresses (void)
00203 {
00204 long i,j,k;
00205
00206 j=adr[0]; adr[0]=0;
00207 for (i=1;i<n+1;i++){
00208 k=adr[i];
00209 adr[i]=j;
00210 j+=k;
00211 }
00212
00213 for (i=0;i<n+1;i++){
00214 adra[i]=adr[i];
00215 }
00216
00217 aux = new long [adr[n]];
00218 memset (aux,0,adr[n]*sizeof(long));
00219 }
00220
00221
00222
00223
00224
00225
00226
00227
00228
00229 void comprow::sort_and_colindex (void)
00230 {
00231 long i,j,k,ii,jj,lj,uj,min,prev;
00232
00233 for (i=0;i<n;i++){
00234 lj=adr[i]; uj=adra[i]; prev=-1;
00235 for (j=lj;j<uj;j++){
00236 min=LONG_MAX;
00237 for (k=j;k<uj;k++){
00238 if (aux[k]<min){
00239 min=aux[k]; ii=k;
00240 }
00241 }
00242 if (min==prev){
00243 uj--; j--;
00244 aux[ii]=aux[uj];
00245 }
00246 else{
00247 jj=aux[j];
00248 aux[j]=min;
00249 aux[ii]=jj;
00250 prev=min;
00251 }
00252 }
00253 adra[i]=uj;
00254 }
00255
00256 j=0;
00257 for (i=0;i<n;i++){
00258 j+=adra[i]-adr[i];
00259 }
00260
00261
00262 negm=j;
00263
00264
00265
00266 printf ("\n\n negm %ld\n\n",negm);
00267
00268
00269
00270
00271
00272 delete [] ci;
00273 ci = new long [negm];
00274 memset (ci,0,negm*sizeof(long));
00275
00276
00277 delete [] a;
00278 a = new double [negm];
00279 memset (a,0,negm*sizeof(double));
00280
00281 ii=0; lj=0; adr[0]=0;
00282 for (i=0;i<n;i++){
00283 uj=adra[i];
00284 for (j=lj;j<uj;j++){
00285 ci[ii]=aux[j]; ii++;
00286 }
00287 lj=adr[i+1];
00288 adr[i+1]=ii;
00289 }
00290
00291 for (i=0;i<=n;i++){
00292 adra[i]=adr[i];
00293 }
00294
00295 delete [] aux;
00296 aux=NULL;
00297 }
00298
00299
00300
00301
00302
00303
00304
00305
00306
00307
00308
00309 void comprow::localize (matrix &b,long *cn)
00310 {
00311 long i,j,k,ii,jj,kk,ll,lk,uk;
00312
00313 for (i=0;i<b.m;i++){
00314 ii=cn[i]-1;
00315 if (ii<0) continue;
00316 ll=0; lk=adr[ii]; uk=adr[ii+1]; kk=lk;
00317 for (j=0;j<b.m;j++){
00318 jj=cn[j]-1;
00319 if (jj<0) continue;
00320 if (ll<jj){
00321 ll=jj;
00322 for (k=kk;k<uk;k++){
00323 if (ci[k]!=jj) continue;
00324 else{
00325 a[k]+=b[i][j];
00326 kk=k;
00327 break;
00328 }
00329 }
00330 }
00331 else{
00332 ll=jj;
00333 for (k=kk;k>=lk;k--){
00334 if (ci[k]!=jj) continue;
00335 else{
00336 a[k]+=b[i][j];
00337 kk=k;
00338 break;
00339 }
00340 }
00341 }
00342 }
00343 }
00344 }
00345
00346
00347
00348
00349
00350
00351
00352
00353
00354
00355
00356 void comprow::localized (double *b,long *cn,long nc)
00357 {
00358 long i,j,k,ii,jj,kk,ll,lk,uk;
00359
00360 for (i=0;i<nc;i++){
00361 ii=cn[i]-1;
00362 if (ii<0) continue;
00363 ll=0; lk=adr[ii]; uk=adr[ii+1]; kk=lk;
00364 for (j=0;j<nc;j++){
00365 jj=cn[j]-1;
00366 if (jj<0) continue;
00367 if (ll<jj){
00368 ll=jj;
00369 for (k=kk;k<uk;k++){
00370 if (ci[k]!=jj) continue;
00371 else{
00372 a[k]+=b[i*nc+j];
00373 kk=k;
00374 break;
00375 }
00376 }
00377 }
00378 else{
00379 ll=jj;
00380 for (k=kk;k>=lk;k--){
00381 if (ci[k]!=jj) continue;
00382 else{
00383 a[k]+=b[i*nc+j];
00384 kk=k;
00385 break;
00386 }
00387 }
00388 }
00389 }
00390 }
00391 }
00392
00393
00394
00395
00396
00397
00398
00399
00400
00401 long comprow::minimize (double limit)
00402 {
00403 long i,j,n1,cor;
00404
00405 cor=-1; n1=0;
00406 for (i=0;i<n;i++){
00407 adr[i]=adra[i]-n1;
00408 for (j=adra[i];j<adra[i+1];j++){
00409 cor++;
00410 if (fabs(a[j])>limit){
00411 a[cor]=a[j];
00412 ci[cor]=ci[j];
00413 }
00414 else{
00415 cor--; n1++;
00416 }
00417 }
00418 }
00419 adr[n]=adra[n]-n1;
00420 negm=adr[n];
00421 delete [] adra;
00422 adra=NULL;
00423 return n1;
00424 }
00425
00426
00427
00428
00429
00430
00431
00432
00433
00434
00435
00436 void comprow::initiate (gtopology *top,long ndof,long mespr)
00437 {
00438 if (status()==NULL){
00439 allocadr (ndof);
00440 numcontr (top);
00441 addresses ();
00442 fillarray (top);
00443 sort_and_colindex ();
00444 }
00445 else{
00446 nullmat ();
00447 }
00448
00449 if (mespr==1) fprintf (stdout,"\n number of matrix entries %ld",negm);
00450 }
00451
00452
00453
00454
00455
00456
00457
00458
00459 void comprow::printmat (FILE *out)
00460 {
00461
00462
00463
00464
00465
00466
00467
00468
00469
00470
00471
00472
00473
00474
00475
00476
00477
00478
00479
00480
00481
00482
00483
00484
00485
00486
00487
00488
00489
00490
00491
00492
00493
00494
00495
00496
00497
00498
00499
00500
00501
00502
00503
00504
00505
00506
00507
00508 long i,j;
00509 FILE *p;
00510
00511 p = fopen ("matice.txt","w");
00512 fprintf (p,"%ld %ld\n",n,negm);
00513
00514 for (i=0;i<n;i++){
00515 for (j=adr[i];j<adr[i+1];j++){
00516
00517
00518 fprintf (p,"%8ld %8ld % 20.15le\n",i,ci[j],a[j]);
00519 }
00520 }
00521 fclose (p);
00522
00523
00524
00525
00526
00527
00528
00529
00530
00531
00532
00533
00534
00535
00536
00537
00538
00539
00540
00541
00542
00543
00544
00545 }
00546
00547
00548
00549
00550
00551
00552
00553
00554
00555 void comprow::printdiag (FILE *out)
00556 {
00557 long i,j;
00558
00559 fprintf (out,"\n\n diagonal entries of matrix stored in the compressed row storage scheme\n");
00560
00561 for (i=0;i<n;i++){
00562 for (j=adr[i];j<adr[i+1];j++){
00563 if (ci[j]==i){
00564 fprintf (out,"%lf\n",a[ci[j]]);
00565 break;
00566 }
00567 }
00568 }
00569 }
00570
00571
00572
00573
00574
00575
00576
00577
00578
00579 double comprow::give_entry (long ri,long rci)
00580 {
00581 long i;
00582 double e=0.0;
00583
00584 for (i=adr[ri];i<adr[ri+1];i++){
00585 if (rci==ci[i]){
00586 e=a[i];
00587 break;
00588 }
00589 }
00590
00591 return e;
00592 }
00593
00594
00595
00596
00597
00598
00599
00600
00601
00602
00603 void comprow::add_entry (double e,long ri,long rci)
00604 {
00605 long i;
00606
00607 for (i=adr[ri];i<adr[ri+1];i++){
00608 if (rci==ci[i]){
00609 a[i]+=e;
00610 break;
00611 }
00612 }
00613 }
00614
00615
00616
00617
00618
00619
00620
00621
00622
00623
00624 void comprow::mxv_cr (double *b,double *c)
00625 {
00626 long i,j,ii,lj,uj;
00627 double s;
00628
00629 for (i=0;i<n;i++){
00630 lj=adr[i]; uj=adr[i+1];
00631 s=0.0;
00632 for (j=lj;j<uj;j++){
00633 ii=ci[j];
00634 s+=a[j]*b[ii];
00635 }
00636 c[i]=s;
00637 }
00638 }
00639
00640
00641
00642
00643
00644
00645
00646
00647
00648
00649
00650
00651 void comprow::mxv_cr_new (double *b,double *c)
00652 {
00653 long i,j,ii,j2,lj,uj,uj2,p1,p2;
00654 double s1,s2;
00655
00656 lj=adr[0];
00657 for (i=0;i<(n-1);i+=2)
00658 {
00659 uj=adr[i+1];
00660 uj2=adr[i+2];
00661 p1=uj-lj;
00662 p2=uj2-uj;
00663 s1=0.0;
00664 s2=0.0;
00665 if (p2>=p1)
00666 {
00667 for (j=lj,j2=uj;j<uj;j++)
00668 {
00669 s1+=a[j]*b[ci[j]];
00670 s2+=a[j2]*b[ci[j2]];
00671 j2++;
00672 }
00673 for (;j2<uj2;j2++)
00674 {
00675 s2+=a[j2]*b[ci[j2]];
00676 }
00677 }
00678 else
00679 {
00680 for (j=lj,j2=uj;j2<uj2;j2++)
00681 {
00682 s1+=a[j]*b[ci[j]];
00683 s2+=a[j2]*b[ci[j2]];
00684 j++;
00685 }
00686 for (;j<uj;j++)
00687 {
00688 s1+=a[j]*b[ci[j]];
00689 }
00690 }
00691 lj=uj2;
00692 c[i]=s1;
00693 c[i+1]=s2;
00694 }
00695 if (n&1)
00696 {
00697 lj=adr[n-1]; uj=adr[n];
00698 s1=0;
00699 for (j=lj;j<uj;j++)
00700 {
00701 ii=ci[j];
00702 s1+=a[j]*b[ii];
00703 }
00704 c[n-1]=s1;
00705 }
00706 }
00707
00708
00709
00710
00711
00712
00713
00714
00715
00716
00717
00718 double comprow::mxv_cr_new2 (double *b,double *c)
00719 {
00720 long i,j,ii,j2,lj,uj,uj2,p1,p2;
00721 double s1,s2,suma;
00722
00723 suma=0.0;
00724 lj=adr[0];
00725 for (i=0;i<(n-1);i+=2)
00726 {
00727 uj=adr[i+1];
00728 uj2=adr[i+2];
00729 p1=uj-lj;
00730 p2=uj2-uj;
00731 s1=0.0;
00732 s2=0.0;
00733 if (p2>=p1)
00734 {
00735 for (j=lj,j2=uj;j<uj;j++)
00736 {
00737 s1+=a[j]*b[ci[j]];
00738 s2+=a[j2]*b[ci[j2]];
00739 j2++;
00740 }
00741 for (;j2<uj2;j2++)
00742 {
00743 s2+=a[j2]*b[ci[j2]];
00744 }
00745 }
00746 else
00747 {
00748 for (j=lj,j2=uj;j2<uj2;j2++)
00749 {
00750 s1+=a[j]*b[ci[j]];
00751 s2+=a[j2]*b[ci[j2]];
00752 j++;
00753 }
00754 for (;j<uj;j++)
00755 {
00756 s1+=a[j]*b[ci[j]];
00757 }
00758 }
00759 lj=uj2;
00760 c[i]=s1;
00761 c[i+1]=s2;
00762 suma+=s1*b[i];
00763 suma+=s2*b[i+1];
00764 }
00765 if (n&1)
00766 {
00767 lj=adr[n-1]; uj=adr[n];
00768 s1=0;
00769 for (j=lj;j<uj;j++)
00770 {
00771 ii=ci[j];
00772 s1+=a[j]*b[ii];
00773 }
00774 c[n-1]=s1;
00775 suma+=s1*b[n-1];
00776 }
00777 return suma;
00778 }
00779
00780
00781
00782
00783
00784
00785
00786
00787
00788
00789 void comprow::mtxv_cr (double *b,double *c)
00790 {
00791 long i,j,ii,lj,uj;
00792 double s;
00793
00794 nullv (c,n);
00795
00796 for (i=0;i<n;i++){
00797 lj=adr[i]; uj=adr[i+1]; s=b[i];
00798 for (j=lj;j<uj;j++){
00799 ii=ci[j];
00800 c[ii]+=a[j]*s;
00801 }
00802 }
00803 }
00804
00805 void comprow::mv_cr15 (double *b,double *c)
00806
00807
00808
00809
00810
00811
00812
00813
00814
00815
00816
00817
00818
00819
00820
00821
00822
00823
00824 {
00825 long i,j,ii,j2,lj,uj,uj2,p1,p2;
00826 double s1,s2;
00827
00828 lj=adr[0];
00829 for (i=0;i<(n-1);i+=2)
00830 {
00831 uj=adr[i+1];
00832 uj2=adr[i+2];
00833 p1=uj-lj;
00834 p2=uj2-uj;
00835 s1=0.0;
00836 s2=0.0;
00837 if (p2>=p1)
00838 {
00839 for (j=lj,j2=uj;j<uj;j++)
00840 {
00841 s1+=a[j]*b[ci[j]];
00842 s2+=a[j2]*b[ci[j2]];
00843 j2++;
00844 }
00845 for (;j2<uj2;j2++)
00846 {
00847 s2+=a[j2]*b[ci[j2]];
00848 }
00849 }
00850 else
00851 {
00852 for (j=lj,j2=uj;j2<uj2;j2++)
00853 {
00854 s1+=a[j]*b[ci[j]];
00855 s2+=a[j2]*b[ci[j2]];
00856 j++;
00857 }
00858 for (;j<uj;j++)
00859 {
00860 s1+=a[j]*b[ci[j]];
00861 }
00862 }
00863 lj=uj2;
00864 c[i]=s1;
00865 c[i+1]=s2;
00866 }
00867 if (n&1)
00868 {
00869 lj=adr[n-1]; uj=adr[n];
00870 s1=0;
00871 for (j=lj;j<uj;j++)
00872 {
00873 ii=ci[j];
00874 s1+=a[j]*b[ii];
00875 }
00876 c[n-1]=s1;
00877 }
00878
00879
00880 }
00881
00882
00883
00884
00885
00886
00887
00888
00889
00890
00891
00892 double comprow::mxv_cr_pom (double *b1,double *c1,double *b2,double *c2)
00893 {
00894 long i,j,ii,lj,uj;
00895 double s1,s2,pom,suma;
00896
00897 suma=0.0;
00898 nullv (c2,n);
00899
00900 for (i=0;i<n;i++)
00901 {
00902 lj=adr[i];
00903 uj=adr[i+1];
00904 s2=b2[i];
00905 s1=0.0;
00906 for (j=lj;j<uj;j++)
00907 {
00908 ii=ci[j];
00909 pom=a[j];
00910 s1+=pom*b1[ii];
00911 c2[ii]+=pom*s2;
00912 }
00913 c1[i]=s1;
00914 suma+=s1*s2;
00915 }
00916 return suma;
00917 }
00918
00919
00920
00921
00922
00923
00924
00925
00926
00927 void comprow::addmat_cr (double c,comprow &cr)
00928 {
00929 long i;
00930 for (i=0;i<negm;i++){
00931 a[i]+=c*cr.a[i];
00932 }
00933 }
00934
00935
00936
00937
00938
00939
00940
00941
00942 void comprow::scalmat_cr (double c)
00943 {
00944 long i;
00945 for (i=0;i<negm;i++){
00946 a[i]*=c;
00947 }
00948 }
00949
00950
00951
00952
00953
00954
00955
00956
00957
00958
00959
00960
00961
00962
00963
00964
00965
00966
00967
00968 void comprow::cg (double *x,double *y,
00969 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_cr (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_cr (d,p);
01009
01010
01011 denom = ss (d,p,n);
01012 if (fabs(denom)<zero){
01013 fprintf (stdout,"\n there is zero denominator (%le) in alpha computation in conjugate gradient method (cr.cpp)\n",fabs(denom));
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 if (i%100==0)
01030 fprintf (stdout,"\n iteration %ld norres %le norres/norrhs %e",i,nom,nom/nory);
01031
01032
01033 if (sqrt(nom/nory)<err) break;
01034
01035
01036
01037 beta = nom/denom;
01038
01039
01040 for (j=0;j<n;j++){
01041 d[j]=beta*d[j]-r[j];
01042 }
01043 }
01044
01045 ani=i; ares=nom;
01046
01047 delete [] p; delete [] r; delete [] d;
01048 }
01049
01050
01051
01052
01053
01054
01055
01056
01057
01058
01059
01060
01061
01062
01063
01064
01065
01066
01067
01068
01069 void comprow::cg_prec (precond &pr,double *x,double *y,
01070 long ni,double err,long &ani,double &ares,double zero,long iv)
01071 {
01072 long i,j;
01073 double nom,denom,nory,alpha,beta;
01074 double *d,*r,*p,*h;
01075
01076 d = new double [n];
01077 r = new double [n];
01078 p = new double [n];
01079 h = new double [n];
01080
01081
01082 if (iv==0){
01083 for (i=0;i<n;i++)
01084 x[i]=0.0;
01085 }
01086
01087
01088 mxv_cr (x,p);
01089
01090
01091 nory=0.0;
01092 for (i=0;i<n;i++){
01093 nory+=y[i]*y[i];
01094 r[i]=y[i]-p[i];
01095 }
01096
01097 if (nory<zero){
01098 print_err("norm of right hand side in conjugate gradient method is smaller than %e", __FILE__, __LINE__, __func__,zero);
01099 ares=nory; ani=0;
01100 return;
01101 }
01102
01103
01104 pr.preconditioning (r,h,r);
01105
01106 nom = ss (r,h,n);
01107
01108
01109 for (i=0;i<n;i++){
01110 d[i]=h[i];
01111 }
01112
01113
01114
01115 for (i=0;i<ni;i++){
01116
01117
01118
01119 mxv_cr (d,p);
01120
01121
01122 denom = ss (d,p,n);
01123 if (fabs(denom)<zero){
01124 fprintf (stdout,"\n there is zero denominator (%le) in alpha computation in conjugate gradient method (cr.cpp)\n",fabs(denom));
01125 break;
01126 }
01127
01128
01129 alpha = nom/denom;
01130
01131
01132 for (j=0;j<n;j++){
01133 x[j]+=alpha*d[j];
01134 r[j]-=alpha*p[j];
01135 }
01136
01137
01138
01139 pr.preconditioning (r,h,r);
01140
01141
01142
01143
01144
01145
01146
01147
01148 denom=nom;
01149 nom = ss (r,h,n);
01150 beta = nom/denom;
01151
01152
01153
01154
01155 if (i%100==0)
01156 fprintf (stdout,"\n iteration %ld norres %le norrhs %le norres/norrhs %e",i,sqrt(nom),sqrt(nory),sqrt(nom/nory));
01157
01158
01159
01160
01161 if (sqrt(nom/nory)<err){
01162 ani=i; ares=nom/nory;
01163 break;
01164 }
01165
01166
01167
01168 for (j=0;j<n;j++){
01169 d[j]=beta*d[j]+h[j];
01170 }
01171 }
01172
01173 ani=i; ares=nom;
01174
01175 delete [] h; delete [] p; delete [] r; delete [] d;
01176 }
01177
01178
01179
01180
01181
01182
01183
01184
01185
01186
01187
01188
01189
01190
01191
01192
01193
01194
01195
01196
01197 void comprow::cg_new (double *x,double *y,
01198 long ni,double err,long &ani,double &ares,double zero,long iv)
01199 {
01200 long i,j;
01201 double nom,denom,nory,alpha,beta,nomx;
01202 double *d,*r,*p;
01203
01204 d = new double [n];
01205 r = new double [n];
01206 p = new double [n];
01207
01208
01209 if (iv==0){
01210 for (i=0;i<n;i++)
01211 x[i]=0.0;
01212 }
01213 mxv_cr_new (x,p);
01214
01215 nory=0.0; nom=0.0;
01216 for (i=0;i<n;i++){
01217 nory+=y[i]*y[i];
01218 r[i]=p[i]-y[i];
01219 nom+=r[i]*r[i];
01220 d[i]=-1.0*r[i];
01221 }
01222
01223 if (nory<zero){
01224 print_err("norm of right hand side in conjugate gradient method is smaller than %e", __FILE__, __LINE__, __func__,zero);
01225 ares=nory; ani=0;
01226 return;
01227 }
01228
01229
01230
01231 for (i=0;i<ni;i++){
01232
01233
01234
01235
01236
01237 denom=mxv_cr_new2(d,p);
01238 if (fabs(denom)<zero){
01239 fprintf (stdout,"\n there is zero denominator in alpha computation in conjugate gradient method (cr.cpp)\n");
01240 break;
01241 }
01242
01243 alpha = nom/denom;
01244
01245
01246 nomx=0.0;
01247 for (j=0;j<n;j++){
01248 r[j]+=alpha*p[j];
01249 x[j]+=alpha*d[j];
01250 nomx+=r[j]*r[j];
01251 }
01252
01253 denom=nom;
01254
01255
01256 nom = nomx;
01257
01258 if (i%100==0)
01259 fprintf (stdout,"\n iteration %ld norres %le norres/norrhs %e",i,nom,nom/nory);
01260
01261
01262 if (nom/nory<err) break;
01263
01264
01265
01266 beta = nom/denom;
01267
01268
01269 for (j=0;j<n;j++){
01270 d[j]=beta*d[j]-r[j];
01271 }
01272 }
01273
01274 ani=i; ares=nom;
01275
01276 delete [] p; delete [] r; delete [] d;
01277 }
01278
01279
01280
01281
01282
01283
01284
01285
01286
01287
01288
01289
01290
01291
01292
01293
01294
01295
01296
01297 void comprow::bicg (double *x,double *y,long ni,double err,long &ani,double &ares,double zero,long iv)
01298 {
01299 long i,j;
01300 double normy,normr,normrr,nom,denom,alpha,beta;
01301 double *r,*rr,*d,*dd,*p,*pp;
01302
01303 r = new double [n]; rr = new double [n];
01304 d = new double [n]; dd = new double [n];
01305 p = new double [n]; pp = new double [n];
01306
01307 normy = ss (y,y,n);
01308
01309 if (normy<zero){
01310 print_err("norm of right hand side in conjugate gradient method is smaller than %e", __FILE__, __LINE__, __func__,zero);
01311 ares=normy; ani=0;
01312 return;
01313 }
01314
01315 if (iv==0){
01316 for (i=0;i<n;i++)
01317 x[i]=0.0;
01318 }
01319
01320 mxv_cr (x,p);
01321 for (i=0;i<n;i++){
01322 r[i]=y[i]-p[i];
01323 d[i]=r[i];
01324 rr[i]=r[i];
01325 dd[i]=rr[i];
01326 }
01327
01328 nom = ss (r,rr,n);
01329
01330 for (i=0;i<ni;i++){
01331 mxv_cr (d,p);
01332 mtxv_cr (dd,pp);
01333
01334 denom = ss (dd,p,n);
01335 alpha = nom/denom;
01336
01337 for (j=0;j<n;j++){
01338 x[j]+=alpha*d[j];
01339 r[j]-=alpha*p[j];
01340 rr[j]-=alpha*pp[j];
01341 }
01342
01343 normr = ss (r,r,n);
01344 normrr = ss (rr,rr,n);
01345 fprintf (stdout,"\n iteration %ld %e %e",i,normr/normy,normrr/normy);
01346 if (normr/normy<err && normrr/normy<err) break;
01347
01348
01349 denom = nom;
01350 nom = ss (r,rr,n);
01351 beta = nom/denom;
01352
01353 for (j=0;j<n;j++){
01354 d[j]=r[j]+beta*d[j];
01355 dd[j]=rr[j]+beta*dd[j];
01356 }
01357 }
01358
01359 ani=i; ares=normr;
01360 delete [] pp; delete [] p;
01361 delete [] dd; delete [] d;
01362 delete [] rr; delete [] r;
01363 }
01364
01365
01366
01367
01368
01369
01370
01371
01372
01373
01374
01375
01376
01377
01378
01379
01380
01381
01382
01383
01384
01385
01386
01387
01388
01389
01390 double comprow::mv_cr15_rev (double *b,double *c)
01391 {
01392 long i,j,j2,lj,uj,uj2,p1,p2;
01393 double s1,s2,soucet=0;
01394
01395 lj=adr[0];
01396 for (i=0;i<n;i+=2)
01397 {
01398 uj=adr[i+1];
01399 uj2=adr[i+2];
01400 p1=uj-lj;
01401 p2=uj2-uj;
01402 s1=0.0;
01403 s2=0.0;
01404 if (p2>=p1)
01405 {
01406 for (j=lj,j2=uj;j<uj;j++)
01407 {
01408 s1+=a[j]*b[ci[j]];
01409 s2+=a[j2]*b[ci[j2]];
01410 j2++;
01411 }
01412 for (;j2<uj2;j2++)
01413 {
01414 s2+=a[j2]*b[ci[j2]];
01415 }
01416 }
01417 else
01418 {
01419 for (j=lj,j2=uj;j2<uj2;j2++)
01420 {
01421 s1+=a[j]*b[ci[j]];
01422 s2+=a[j2]*b[ci[j2]];
01423 j++;
01424 }
01425 for (;j<uj;j++)
01426 {
01427 s1+=a[j]*b[ci[j]];
01428 }
01429 }
01430 lj=uj2;
01431 c[i]=s1;
01432 c[i+1]=s2;
01433 soucet+=s1*b[i]+s2*b[i+1];
01434 }
01435 if(n&1)
01436 {
01437 s1=0.0;
01438 for (j=adr[n-1];j<adr[n];j++)
01439 {
01440 s1+=a[j]*b[ci[j]];
01441 }
01442 c[n-1]=s1;
01443
01444 soucet+=s1*b[n-1];
01445 }
01446 return soucet;
01447 }
01448
01449
01450
01451
01452
01453
01454
01455
01456
01457
01458
01459
01460
01461
01462
01463
01464
01465
01466
01467
01468
01469 void comprow::cg_cr_rev (double *x,double *y,
01470 long ni,double err,long &ani,double &ares,double limit,long iv)
01471 {
01472 long i,j;
01473 double nom,denom,nory,alpha,beta,ener;
01474 double *d,*r,*p;
01475
01476 d = new double [n];
01477 r = new double [n];
01478 p = new double [n];
01479
01480
01481 if (iv==0){
01482 for (i=0;i<n;i++)
01483 x[i]=0.0;
01484 }
01485 mv_cr15_rev (x,p);
01486
01487 nory=0.0; nom=0.0;
01488 for (i=0;i<n;i++){
01489 nory+=y[i]*y[i];
01490 r[i]=0.0-y[i];
01491 d[i]=y[i];
01492 nom+=r[i]*r[i];
01493 }
01494
01495 if (nory<limit){
01496 printf ("\n\n norma vektoru prave strany v metode sdruzenych gradientu");
01497 printf ("\n je mensi nez %e",limit);
01498 ares=nory; ani=0;
01499 return;
01500 }
01501
01502 for (i=0;i<ni;i++)
01503 {
01504
01505 denom=mv_cr15_rev (d,p);
01506
01507
01508 if (fabs(denom)<limit){
01509 printf ("\n V metode sdruzenych gradientu je nulovy jmenovatel u alfa");
01510 abort ();
01511 }
01512
01513 alpha = nom/denom;
01514 ener=0;
01515
01516 for (j=0;j<n;j++){
01517 x[j]+=alpha*d[j];
01518 r[j]+=alpha*p[j];
01519 ener+=r[j]*r[j];
01520 }
01521
01522
01523
01524
01525
01526
01527 denom=nom;
01528
01529
01530 nom=ener;
01531
01532
01533
01534 if (nom/nory<err){
01535 printf ("\n iterace %ld %e",i,nom/nory);
01536
01537 break;
01538 }
01539 if (fabs(nom)<limit){
01540 printf ("\n iterace %ld %e",i,nom/nory);
01541
01542 break;
01543 }
01544
01545 beta = nom/denom;
01546
01547
01548 for (j=n-1;j>=0;j--){
01549 d[j]=beta*d[j]-r[j];
01550 }
01551 }
01552
01553 ani=i; ares=nom;
01554
01555 delete [] p; delete [] r; delete [] d;
01556 }
01557
01558
01559
01560
01561
01562
01563
01564
01565
01566
01567
01568
01569
01570
01571
01572
01573
01574
01575
01576
01577 void comprow::bicg_new (double *x,double *y,long ni,double err,long &ani,double &ares,double zero,long iv)
01578 {
01579 long i,j;
01580 double normy,normr,normrr,nom,denom,alpha,beta,nomx;
01581 double *r,*rr,*d,*dd,*p,*pp;
01582
01583 r = new double [n]; rr = new double [n];
01584 d = new double [n]; dd = new double [n];
01585 p = new double [n]; pp = new double [n];
01586
01587 normy = ss (y,y,n);
01588
01589 if (normy<zero){
01590 print_err("norm of right hand side in conjugate gradient method is smaller than %e", __FILE__, __LINE__, __func__,zero);
01591 ares=normy; ani=0;
01592 return;
01593 }
01594
01595 if (iv==0){
01596 for (i=0;i<n;i++)
01597 x[i]=0.0;
01598 }
01599
01600 mxv_cr (x,p);
01601 for (i=0;i<n;i++){
01602 r[i]=y[i]-p[i];
01603 d[i]=r[i];
01604 rr[i]=r[i];
01605 dd[i]=rr[i];
01606 }
01607
01608 nom = ss (r,rr,n);
01609
01610 for (i=0;i<ni;i++){
01611
01612
01613
01614 denom=mxv_cr_pom (d,p,dd,pp);
01615
01616 alpha = nom/denom;
01617 normr = normrr = nomx = 0.0;
01618 for (j=0;j<n;j++){
01619 x[j]+=alpha*d[j];
01620 r[j]-=alpha*p[j];
01621 normr += r[j]*r[j];
01622 rr[j]-=alpha*pp[j];
01623 normrr+= rr[j]*rr[j];
01624 nomx+=r[j]*rr[j];
01625 }
01626
01627
01628
01629 fprintf (stdout,"\n iteration %ld %e %e",i,normr/normy,normrr/normy);
01630 if (normr/normy<err && normrr/normy<err) break;
01631
01632
01633 denom = nom;
01634
01635 nom=nomx;
01636 beta = nom/denom;
01637
01638 for (j=0;j<n;j++){
01639 d[j]=r[j]+beta*d[j];
01640 dd[j]=rr[j]+beta*dd[j];
01641 }
01642 }
01643
01644 ani=i; ares=normr;
01645 delete [] pp; delete [] p;
01646 delete [] dd; delete [] d;
01647 delete [] rr; delete [] r;
01648 }
01649
01650
01651
01652
01653
01654
01655
01656
01657 void comprow::copy_cr (comprow &cr)
01658 {
01659 long i;
01660
01661 cr.n=n;
01662 cr.negm=negm;
01663
01664 if (cr.adr!=NULL){
01665 delete [] cr.adr;
01666 }
01667 cr.adr=new long [n+1];
01668
01669 if (cr.ci!=NULL){
01670 delete [] cr.ci;
01671 }
01672 cr.ci=new long [negm];
01673
01674 if (cr.a!=NULL){
01675 delete [] cr.a;
01676 }
01677 cr.a=new double [negm];
01678
01679 for (i=0;i<=n;i++){
01680 cr.adr[i]=adr[i];
01681 }
01682 for (i=0;i<negm;i++){
01683 cr.a[i]=a[i];
01684 cr.ci[i]=ci[i];
01685 }
01686 }
01687
01688
01689
01690
01691
01692
01693
01694
01695
01696
01697
01698
01699 void comprow::fill_data (long nsub,long *smadr,long *smci,double *sma)
01700 {
01701 long i;
01702
01703
01704 n=nsub;
01705
01706
01707 if (adr!=NULL){
01708 delete [] adr;
01709 }
01710 adr = new long [n+1];
01711
01712 for (i=0;i<n+1;i++){
01713 adr[i]=smadr[i];
01714 }
01715
01716
01717 negm = adr[n];
01718
01719
01720 if (ci!=NULL){
01721 delete [] ci;
01722 }
01723 ci = new long [negm];
01724
01725 for (i=0;i<negm;i++){
01726 ci[i]=smci[i];
01727 }
01728
01729
01730 if (a!=NULL){
01731 delete [] a;
01732 }
01733 a = new double [negm];
01734
01735 for (i=0;i<negm;i++){
01736 a[i]=sma[i];
01737 }
01738
01739
01740
01741 }
01742
01743
01744
01745
01746
01747
01748
01749
01750
01751
01752
01753 void comprow::select_submatrix (long *li,long nsub,comprow *smcr)
01754 {
01755 long i,j,k,m,ii,minm,maxm,mins,maxs;
01756 long *smadr,*smci,*newnumb;
01757 double *sma;
01758
01759 smadr = new long [nsub+1];
01760
01761
01762 mins=li[0]-1;
01763
01764 maxs=li[nsub-1]-1;
01765
01766
01767
01768
01769 for (i=0;i<nsub;i++){
01770 smadr[i]=0;
01771
01772 m=li[i]-1;
01773
01774 minm=ci[adr[m]];
01775
01776 maxm=ci[adr[m+1]-1];
01777
01778 if (maxm<mins)
01779 continue;
01780 if (maxs<minm)
01781 continue;
01782
01783 j=adr[m]; k=0;
01784 do{
01785 if (ci[j]==li[k]-1){
01786 smadr[i]++;
01787 j++; k++;
01788 }
01789 else{
01790 if (ci[j]<li[k]-1){
01791 j++;
01792 }
01793 else{
01794 k++;
01795 }
01796 }
01797 }while (j<adr[m+1] && k<nsub);
01798 }
01799
01800
01801
01802 j=0;
01803 for (i=0;i<nsub;i++){
01804 j+=smadr[i];
01805 }
01806
01807 smci = new long [j];
01808 sma = new double [j];
01809
01810 ii=0;
01811 for (i=0;i<nsub;i++){
01812
01813 m=li[i]-1;
01814
01815 minm=ci[adr[m]];
01816
01817 maxm=ci[adr[m+1]-1];
01818
01819 if (maxm<mins)
01820 continue;
01821 if (maxs<minm)
01822 continue;
01823
01824 j=adr[m]; k=0;
01825 do{
01826 if (ci[j]==li[k]-1){
01827 sma[ii]=a[j];
01828 smci[ii]=ci[j];
01829 ii++;
01830 j++; k++;
01831 }
01832 else{
01833 if (ci[j]<li[k]-1){
01834 j++;
01835 }
01836 else{
01837 k++;
01838 }
01839 }
01840 }while (j<adr[m+1] && k<nsub);
01841
01842 }
01843
01844 m=smadr[0]; smadr[0]=0;
01845 for (j=1;j<nsub+1;j++){
01846 k=smadr[j];
01847 smadr[j]=m;
01848 m+=k;
01849 }
01850
01851
01852
01853
01854
01855
01856
01857
01858 newnumb = new long [maxs+1];
01859 for (i=0;i<maxs+1;i++){
01860 newnumb[i]=-1;
01861 }
01862 for (i=0;i<nsub;i++){
01863 j=li[i]-1;
01864 newnumb[j]=1;
01865 }
01866 j=0;
01867 for (i=0;i<maxs+1;i++){
01868 if (newnumb[i]==1){
01869 newnumb[i]=j;
01870 j++;
01871 }
01872 }
01873
01874 for (i=0;i<smadr[nsub];i++){
01875 j=smci[i];
01876 smci[i]=newnumb[j];
01877 }
01878
01879
01880
01881 smcr->fill_data (nsub,smadr,smci,sma);
01882
01883 delete [] sma;
01884 delete [] smci;
01885 delete [] smadr;
01886 delete [] aux;
01887 aux = NULL;
01888 }
01889
01890
01891
01892
01893
01894
01895
01896
01897
01898
01899
01900
01901 void comprow::select_submatrix (long *li,long nsub,densemat *smdm)
01902 {
01903 long i,j,k,ii,jj,kk;
01904
01905 smdm->n = nsub;
01906 smdm->a = new double [nsub*nsub];
01907 for (i=0;i<nsub;i++){
01908 for (j=0;j<nsub;j++){
01909 smdm->a[i*nsub+j]=0.0;
01910 }
01911 }
01912
01913 for (i=0;i<nsub;i++){
01914
01915 ii=li[i]-1;
01916 for (j=adr[ii];j<adr[ii+1];j++){
01917
01918 jj=ci[j]; kk=-1;
01919 for (k=0;k<nsub;k++){
01920 if (li[k]-1==jj){
01921 kk=k;
01922 break;
01923 }
01924 }
01925 if (kk>-1)
01926 smdm->a[i*nsub+k]=a[j];
01927 }
01928 }
01929
01930 }
01931
01932
01933
01934
01935
01936
01937
01938
01939
01940
01941 void comprow::crxcv_cv (compvect *cvi,compvect *cvo)
01942 {
01943 long i,j,k;
01944 long minm,maxm,minv,maxv,ncm,ncv;
01945 long *ind;
01946 double s,*b,*c;
01947
01948 if (cvi->ordering == 0){
01949 print_err("compressed sparse vector is not ordered increasingly", __FILE__, __LINE__, __func__);
01950 }
01951
01952
01953 c = new double [n];
01954
01955
01956 ncv=cvi->nz;
01957
01958 minv=cvi->ind[0];
01959
01960 maxv=cvi->ind[ncv-1];
01961
01962 ind=cvi->ind;
01963 b=cvi->a;
01964
01965 for (i=0;i<n;i++){
01966 c[i]=0.0;
01967
01968
01969 ncm=adr[i+1]-adr[i];
01970
01971 minm=ci[adr[i]];
01972
01973 maxm=ci[adr[i+1]-1];
01974
01975 if (maxm<minv)
01976 continue;
01977 if (maxv<minm)
01978 continue;
01979
01980 j=adr[i]; k=0;
01981 s=0.0;
01982 do{
01983 if (ci[j]==ind[k]){
01984 s+=a[j]*b[k];
01985 j++; k++;
01986 }
01987 else{
01988 if (ci[j]<ind[k]){
01989 j++;
01990 }
01991 else{
01992 k++;
01993 }
01994 }
01995 }while (j<adr[i+1] && k<ncv);
01996 c[i]=s;
01997 }
01998
01999
02000 cvo->setup_vector (c,n);
02001
02002 delete [] c;
02003 }
02004
02005
02006
02007
02008
02009
02010 long comprow::give_negm ()
02011 {
02012 return negm;
02013 }
02014
02015
02016
02017
02018
02019
02020
02021
02022
02023
02024 double comprow::estim_spect_radius ()
02025 {
02026 long i,j;
02027 double s,sr;
02028
02029
02030 sr=0.0;
02031 for (i=0;i<n;i++){
02032 s=0.0;
02033 for (j=adr[i];j<adr[i+1];j++){
02034 s+=fabs(a[j]);
02035 }
02036 if (s>sr)
02037 sr=s;
02038 }
02039
02040 return sr;
02041 }