00001 #include "scr.h"
00002 #include <limits.h>
00003 #include <math.h>
00004 #include <string.h>
00005
00006
00007 symcomprow::symcomprow (void)
00008 {
00009 limit = 0.0;
00010 adra=NULL;
00011 ci=NULL;
00012 a=NULL;
00013 aux=NULL;
00014 incdec=NULL;
00015 decompid=0;
00016 adr = NULL;
00017 }
00018
00019
00020
00021 symcomprow::~symcomprow (void)
00022 {
00023 delete [] adr;
00024 delete [] ci;
00025 delete [] a;
00026 delete [] incdec;
00027 delete [] adra;
00028 delete [] aux;
00029 }
00030
00031
00032
00033
00034
00035
00036
00037 long symcomprow::decomp ()
00038 {
00039 return decompid;
00040 }
00041
00042
00043
00044
00045
00046
00047 void symcomprow::changedecomp ()
00048 {
00049 if (decompid==0) decompid=1;
00050 else decompid=0;
00051 }
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061 void symcomprow::allocadr (long m)
00062 {
00063 n=m;
00064
00065 adr = new long [n+1];
00066 adra = new long [n+1];
00067
00068 memset (adr,0,(n+1)*sizeof(long));
00069 memset (adra,0,(n+1)*sizeof(long));
00070 }
00071
00072
00073
00074
00075
00076
00077
00078 double* symcomprow::status ()
00079 {
00080 return a;
00081 }
00082
00083
00084
00085
00086
00087
00088 void symcomprow::nullmat ()
00089 {
00090 long i;
00091
00092 for (i=0;i<negm;i++){
00093 a[i]=0.0;
00094 }
00095 }
00096
00097
00098
00099
00100
00101
00102
00103
00104
00105
00106 void symcomprow::numcontr_elem (long *cn,long ndofe)
00107 {
00108 long i,j,ii,jj;
00109
00110 for (i=0;i<ndofe;i++){
00111 ii=cn[i]-1;
00112 if (ii<0) continue;
00113 for (j=0;j<ndofe;j++){
00114 jj=cn[j]-1;
00115 if (jj<0) continue;
00116 if (ii<jj) continue;
00117 adr[ii]++;
00118 }
00119 }
00120 }
00121
00122
00123
00124
00125
00126
00127
00128
00129
00130 void symcomprow::numcontr (gtopology *top)
00131 {
00132 long i,ne,ndofe,*cn;
00133 ne=top->ne;
00134 for (i=0;i<ne;i++){
00135 ndofe=top->give_ndofe (i);
00136 cn = new long [ndofe];
00137 top->give_code_numbers (i,cn);
00138 numcontr_elem (cn,ndofe);
00139 delete [] cn;
00140 }
00141 }
00142
00143
00144
00145
00146
00147
00148
00149
00150
00151
00152
00153 void symcomprow::fillarray_elem (long *cn,long ndofe)
00154 {
00155 long i,j,ii,jj;
00156
00157 for (i=0;i<ndofe;i++){
00158 ii=cn[i]-1;
00159 if (ii<0) continue;
00160 for (j=0;j<ndofe;j++){
00161 jj=cn[j]-1;
00162 if (jj<0) continue;
00163 if (ii<jj) continue;
00164 aux[adra[ii]]=jj;
00165 adra[ii]++;
00166 }
00167 }
00168 }
00169
00170
00171
00172
00173
00174
00175
00176
00177
00178
00179 void symcomprow::fillarray (gtopology *top)
00180 {
00181 long i,ne,ndofe,*cn;
00182 ne=top->ne;
00183 for (i=0;i<ne;i++){
00184 ndofe=top->give_ndofe (i);
00185 cn = new long [ndofe];
00186 top->give_code_numbers (i,cn);
00187 fillarray_elem (cn,ndofe);
00188 delete [] cn;
00189 }
00190 }
00191
00192
00193
00194
00195
00196
00197 void symcomprow::addresses (void)
00198 {
00199 long i,j,k;
00200
00201 j=adr[0]; adr[0]=0;
00202 for (i=1;i<n+1;i++){
00203 k=adr[i];
00204 adr[i]=j;
00205 j+=k;
00206 }
00207
00208 for (i=0;i<n+1;i++){
00209 adra[i]=adr[i];
00210 }
00211
00212 aux = new long [adr[n]];
00213 memset (aux,0,adr[n]*sizeof(long));
00214 }
00215
00216
00217
00218
00219
00220
00221
00222
00223
00224 void symcomprow::sort_and_colindex (void)
00225 {
00226 long i,j,k,ii,jj,lj,uj,min,prev;
00227
00228
00229
00230
00231
00232 for (i=0;i<n;i++){
00233 lj=adr[i]; uj=adra[i]; prev=-1;
00234 for (j=lj;j<uj;j++){
00235 min=LONG_MAX;
00236 for (k=j;k<uj;k++){
00237 if (aux[k]<min){
00238 min=aux[k]; ii=k;
00239 }
00240 }
00241 if (min==prev){
00242 uj--; j--;
00243 aux[ii]=aux[uj];
00244 }
00245 else{
00246 jj=aux[j];
00247 aux[j]=min;
00248 aux[ii]=jj;
00249 prev=min;
00250 }
00251 }
00252 adra[i]=uj;
00253 }
00254
00255 j=0;
00256 for (i=0;i<n;i++){
00257 j+=adra[i]-adr[i];
00258 }
00259
00260
00261 negm=j;
00262
00263
00264 delete [] ci;
00265
00266 ci = new long [negm];
00267 memset (ci,0,negm*sizeof(long));
00268
00269
00270 delete [] a;
00271 a = new double [negm];
00272 memset (a,0,negm*sizeof(double));
00273
00274
00275 ii=0; lj=0; adr[0]=0;
00276 for (i=0;i<n;i++){
00277 uj=adra[i];
00278 for (j=lj;j<uj;j++){
00279 ci[ii]=aux[j]; ii++;
00280 }
00281 lj=adr[i+1];
00282 adr[i+1]=ii;
00283 }
00284
00285 for (i=0;i<=n;i++){
00286 adra[i]=adr[i];
00287 }
00288
00289
00290 delete [] aux;
00291 aux = NULL;
00292
00293 }
00294
00295
00296
00297
00298
00299
00300
00301
00302
00303 void symcomprow::sort_and_colindex_tko (void)
00304 {
00305 long i,j,ii,lj,uj,nd;
00306
00307
00308
00309
00310
00311 long *udof = new long[n];
00312
00313 for (i=0;i<n;i++)
00314 {
00315 memset(udof, 0, sizeof(*udof)*n);
00316
00317 for (j=adr[i]; j<adra[i]; j++)
00318 udof[aux[j]] = 1;
00319
00320 nd=0;
00321 lj = adr[i];
00322 for (j=0; j<n; j++)
00323 {
00324 if (udof[j])
00325 {
00326 aux[lj+nd] = j;
00327 nd++;
00328 }
00329 }
00330 adra[i] = adr[i]+nd;
00331 }
00332
00333 delete [] udof;
00334
00335 j=0;
00336 for (i=0;i<n;i++){
00337 j+=adra[i]-adr[i];
00338 }
00339
00340
00341 negm=j;
00342
00343
00344 delete [] ci;
00345
00346 ci = new long [negm];
00347 memset (ci,0,negm*sizeof(long));
00348
00349
00350 delete [] a;
00351 a = new double [negm];
00352 memset (a,0,negm*sizeof(double));
00353
00354
00355 ii=0; lj=0; adr[0]=0;
00356 for (i=0;i<n;i++){
00357 uj=adra[i];
00358 for (j=lj;j<uj;j++){
00359 ci[ii]=aux[j]; ii++;
00360 }
00361 lj=adr[i+1];
00362 adr[i+1]=ii;
00363 }
00364
00365 for (i=0;i<=n;i++){
00366 adra[i]=adr[i];
00367 }
00368
00369
00370 delete [] aux;
00371 aux = NULL;
00372
00373 }
00374
00375
00376
00377
00378
00379
00380
00381
00382
00383
00384
00385
00386 void symcomprow::localize (matrix &b,long *cn)
00387 {
00388 long i,j,k,ii,jj,kk,lk,uk,ll;
00389
00390 for (i=0;i<b.m;i++){
00391 ii=cn[i]-1;
00392 if (ii<0) continue;
00393 ll=0; lk=adr[ii]; uk=adr[ii+1]; kk=lk;
00394 for (j=0;j<b.m;j++){
00395 jj=cn[j]-1;
00396 if (jj<0) continue;
00397 if (ii<jj) continue;
00398 if (ll<jj){
00399 ll=jj;
00400 for (k=kk;k<uk;k++){
00401 if (ci[k]!=jj) continue;
00402 else{
00403 a[k]+=b[i][j];
00404 kk=k;
00405 break;
00406 }
00407 }
00408 }
00409 else{
00410 ll=jj;
00411 for (k=kk;k>=lk;k--){
00412 if (ci[k]!=jj) continue;
00413 else{
00414 a[k]+=b[i][j];
00415 kk=k;
00416 break;
00417 }
00418 }
00419 }
00420 }
00421 }
00422 }
00423
00424
00425
00426
00427
00428
00429
00430
00431
00432
00433
00434 void symcomprow::localized (double *b,long *cn,long nc)
00435 {
00436 long i,j,k,ii,jj,kk,lk,uk,ll;
00437
00438 for (i=0;i<nc;i++){
00439 ii=cn[i]-1;
00440 if (ii<0) continue;
00441 ll=0; lk=adr[ii]; uk=adr[ii+1]; kk=lk;
00442 for (j=0;j<nc;j++){
00443 jj=cn[j]-1;
00444 if (jj<0) continue;
00445 if (ii<jj) continue;
00446 if (ll<jj){
00447 ll=jj;
00448 for (k=kk;k<uk;k++){
00449 if (ci[k]!=jj) continue;
00450 else{
00451 a[k]+=b[i*nc+j];
00452 kk=k;
00453 break;
00454 }
00455 }
00456 }
00457 else{
00458 ll=jj;
00459 for (k=kk;k>=lk;k--){
00460 if (ci[k]!=jj) continue;
00461 else{
00462 a[k]+=b[i*nc+j];
00463 kk=k;
00464 break;
00465 }
00466 }
00467 }
00468 }
00469 }
00470 }
00471
00472
00473
00474
00475
00476
00477
00478
00479
00480
00481 long symcomprow::minimize (double limit)
00482 {
00483 long i,j,n1,cor;
00484
00485 cor=-1; n1=0;
00486 for (i=0;i<n;i++){
00487 adr[i]=adra[i]-n1;
00488 for (j=adra[i];j<adra[i+1];j++){
00489 cor++;
00490 if (fabs(a[j])>limit){
00491 a[cor]=a[j];
00492 ci[cor]=ci[j];
00493 }
00494 else{
00495 cor--; n1++;
00496 }
00497 }
00498 }
00499 adr[n]=adra[n]-n1;
00500 negm=adr[n];
00501 delete [] adra;
00502 adra = NULL;
00503 return n1;
00504 }
00505
00506
00507
00508
00509
00510
00511
00512
00513
00514
00515
00516 void symcomprow::initiate (gtopology *top,long ndof,long mespr)
00517 {
00518 if (status()==NULL){
00519 allocadr (ndof);
00520 numcontr (top);
00521 addresses ();
00522 fillarray (top);
00523 sort_and_colindex_tko ();
00524 }
00525 else{
00526 nullmat ();
00527 }
00528
00529 if (mespr==1) fprintf (stdout,"\n number of matrix entries (kontrola) %ld",negm);
00530 }
00531
00532
00533
00534
00535
00536
00537
00538
00539 double symcomprow::give_entry (long ri,long rci)
00540 {
00541 long i;
00542 double e=0.0;
00543
00544 if (ri>rci){
00545 for (i=adr[ri];i<adr[ri+1];i++){
00546 if (rci==ci[i]){
00547 e=a[i];
00548 break;
00549 }
00550 }
00551 }
00552 else{
00553 i=ri;
00554 ri=rci;
00555 rci=i;
00556 for (i=adr[ri];i<adr[ri+1];i++){
00557 if (rci==ci[i]){
00558 e=a[i];
00559 break;
00560 }
00561 }
00562 }
00563
00564 return e;
00565 }
00566
00567
00568
00569
00570
00571
00572
00573
00574
00575 void symcomprow::add_entry (double e,long ri,long rci)
00576 {
00577 long i;
00578
00579 if (ri>rci){
00580 for (i=adr[ri];i<adr[ri+1];i++){
00581 if (rci==ci[i]){
00582 a[i]+=e;
00583 break;
00584 }
00585 }
00586 }
00587 else{
00588 i=ri;
00589 ri=rci;
00590 rci=i;
00591 for (i=adr[ri];i<adr[ri+1];i++){
00592 if (rci==ci[i]){
00593 a[i]+=e;
00594 break;
00595 }
00596 }
00597 }
00598
00599 }
00600
00601
00602
00603
00604
00605
00606
00607
00608
00609
00610 void symcomprow::mxv_scr (double *b,double *c)
00611 {
00612 long i,j,ii,lj,uj;
00613 double s,d;
00614
00615 for (i=0;i<n;i++){
00616 lj=adr[i]; uj=adr[i+1];
00617 s=0.0; d=b[i];
00618 for (j=lj;j<uj;j++){
00619 ii=ci[j];
00620 s+=a[j]*b[ii];
00621 c[ii]+=a[j]*d;
00622 }
00623 c[i]=s;
00624 }
00625 }
00626
00627
00628
00629
00630
00631
00632
00633
00634
00635 void symcomprow::addmat_scr (double c,symcomprow &scr)
00636 {
00637 long i;
00638 for (i=0;i<negm;i++){
00639 a[i]+=c*scr.a[i];
00640 }
00641 }
00642
00643
00644
00645
00646
00647
00648
00649
00650 void symcomprow::scalmat_scr (double c)
00651 {
00652 long i;
00653 for (i=0;i<negm;i++){
00654 a[i]*=c;
00655 }
00656 }
00657
00658
00659
00660
00661
00662
00663
00664
00665
00666
00667
00668
00669
00670
00671
00672
00673
00674
00675
00676 void symcomprow::cg (double *x,double *y,long ni,double err,long &ani,double &ares,double zero,long iv)
00677 {
00678 long i,j;
00679 double nom,denom,nory,alpha,beta;
00680 double *d,*r,*p;
00681
00682 fprintf (stdout,"NDOF %ld\n",n);
00683
00684 d = new double [n];
00685 r = new double [n];
00686 p = new double [n];
00687
00688
00689 if (iv==0){
00690 for (i=0;i<n;i++)
00691 x[i]=0.0;
00692 }
00693
00694 mxv_scr (x,p);
00695
00696
00697
00698
00699
00700 nory=0.0; nom=0.0;
00701 for (i=0;i<n;i++){
00702 nory+=y[i]*y[i];
00703 r[i]=y[i]-p[i];
00704 nom+=r[i]*r[i];
00705 d[i]=r[i];
00706 }
00707
00708 fprintf (stdout," SCR nory %le\n",nory);
00709 fprintf (stdout," SCR nom %le\n",nom);
00710
00711 if (nory<zero){
00712 print_warning("norm of right hand side in conjugate gradient method is smaller than %e", __FILE__, __LINE__, __func__,zero);
00713 ares=nory; ani=0;
00714 for (i=0;i<n;i++)
00715 x[i]=0.0;
00716
00717 delete [] p; delete [] r; delete [] d;
00718
00719 return;
00720 }
00721
00722
00723
00724 for (i=0;i<ni;i++){
00725
00726
00727 mxv_scr (d,p);
00728
00729 denom = ss (d,p,n);
00730 if (fabs(denom)<zero){
00731 print_warning("there is zero denominator in alpha computation in conjugate gradient method", __FILE__, __LINE__, __func__);
00732 break;
00733 }
00734
00735 alpha = nom/denom;
00736
00737
00738 for (j=0;j<n;j++){
00739 x[j]+=alpha*d[j];
00740 r[j]-=alpha*p[j];
00741 }
00742
00743 denom=nom;
00744
00745 nom = ss (r,r,n);
00746
00747
00748
00749 fprintf (stdout,"\n iteration %ld norres/norrhs %e",i,nom/nory);
00750
00751
00752 if (nom/nory<err){
00753 fprintf (stdout,"\n\n Last iteration step %ld norres/norrhs %e\n\n",i,nom/nory);
00754 break;
00755 }
00756
00757 if (fabs(nom)<zero){
00758 fprintf (stdout,"\n\n Last iteration step %ld norres/norrhs %e\n\n",i,nom/nory);
00759 break;
00760 }
00761
00762 beta = nom/denom;
00763
00764
00765 for (j=0;j<n;j++){
00766 d[j]=beta*d[j]+r[j];
00767 }
00768 }
00769
00770 ani=i; ares=nom;
00771
00772 delete [] p; delete [] r; delete [] d;
00773 }
00774
00775
00776
00777
00778
00779
00780
00781
00782
00783
00784
00785
00786
00787 void symcomprow::prec_diag_scr (double *x,double *y)
00788 {
00789 long i;
00790
00791 for (i=0;i<n;i++){
00792 x[i]=y[i]/a[adr[i+1]-1];
00793 }
00794 }
00795
00796
00797
00798
00799
00800
00801
00802
00803
00804
00805
00806
00807
00808
00809 void symcomprow::prec_ssor_scr (double *x,double *y,double omega)
00810 {
00811 long i,j,uj;
00812 double s,*v;
00813
00814 v = new double [n];
00815 memset (v,0,n*sizeof(double));
00816
00817 for (i=0;i<n;i++){
00818 uj=adr[i+1]-1; s=0.0;
00819 for (j=adr[i];j<uj;j++){
00820 s+=a[j]*v[ci[j]];
00821 }
00822 v[i]=(y[i]-s)/a[uj]*omega;
00823 }
00824
00825 for (i=0;i<n;i++){
00826 v[i]=v[i]*a[adr[i+1]-1]/omega*(2.0-omega);
00827 }
00828
00829 for (i=n-1;i>-1;i--){
00830 x[i]=v[i]*omega/a[adr[i+1]-1];
00831 uj=adr[i+1]-1; s=x[i];
00832 for (j=adr[i];j<uj;j++){
00833 v[ci[j]]-=a[j]*s;
00834 }
00835 }
00836
00837 delete [] v;
00838 }
00839
00840
00841
00842
00843
00844
00845
00846
00847
00848
00849 double symcomprow::rows_mult_ldl_scr (long m,long o)
00850 {
00851 long i,j,ii,jj;
00852 double s;
00853
00854 if (m<o){
00855 print_err("wrong row numbers", __FILE__, __LINE__, __func__);
00856 }
00857
00858 s=0.0; ii=adr[m]; jj=adr[o];
00859 i=ci[ii]; j=ci[jj];
00860 while (j<o){
00861 if (i==j) { s+=incdec[ii]*incdec[jj]*incdec[adr[j+1]-1]; ii++; jj++; }
00862 if (i<j) { ii++; }
00863 if (i>j) { jj++; }
00864 i=ci[ii]; j=ci[jj];
00865 }
00866
00867 return s;
00868 }
00869
00870
00871
00872
00873
00874
00875
00876
00877
00878
00879
00880 void symcomprow::incomplete_ldl (double gamma)
00881 {
00882 long i,j,ii,lj,uj;
00883 double s;
00884
00885
00886 incdec = new double [negm];
00887
00888 for (i=0;i<n;i++){
00889
00890 ii=adr[i]; lj=ci[ii]; uj=adr[i+1]-1; incdec[uj]=a[uj];
00891 for (j=lj;j<i;j++){
00892 s = rows_mult_ldl_scr (i,j);
00893 if (j==ci[ii]) { incdec[ii]=(a[ii]-s)/incdec[adr[j+1]-1]; ii++; }
00894 else { incdec[uj]-=gamma*s; }
00895 }
00896
00897
00898 lj=adr[i]; s=0.0;
00899 for (j=lj;j<uj;j++){
00900 s+=incdec[j]*incdec[j]*incdec[adr[ci[j]+1]-1];
00901 }
00902 incdec[uj]-=s;
00903 }
00904
00905 }
00906
00907
00908
00909
00910
00911
00912
00913
00914
00915
00916
00917
00918
00919 void symcomprow::prec_ildl_scr (double *x,double *y)
00920 {
00921 long i,j,lj,uj;
00922 double s,*p;
00923
00924 p = new double [n];
00925
00926
00927 for (i=0;i<n;i++){
00928 s=0.0; lj=adr[i]; uj=adr[i+1]-1;
00929 for (j=lj;j<uj;j++){
00930 s+=incdec[j]*p[ci[j]];
00931 }
00932 p[i]=y[i]-s;
00933 }
00934
00935
00936 for (i=0;i<n;i++){
00937 p[i]/=incdec[adr[i+1]-1];
00938 }
00939
00940
00941 for (i=n-1;i>-1;i--){
00942 lj=adr[i]; uj=adr[i+1]-1;
00943 x[i]=p[i]; s=x[i];
00944 for (j=lj;j<uj;j++){
00945 p[ci[j]]-=incdec[j]*s;
00946 }
00947 }
00948
00949 delete [] p;
00950 }
00951
00952
00953
00954
00955
00956
00957
00958
00959
00960
00961
00962
00963
00964
00965
00966
00967
00968
00969
00970
00971
00972
00973
00974
00975
00976
00977
00978
00979
00980
00981
00982
00983 void symcomprow::cg_prec (double *x,double *y,long ni,double err,long &ani,double &ares,
00984 double zero,long iv,long tprec,double par,ISolver *sdirect)
00985 {
00986 long i,j;
00987 double nom,denom,nory,alpha,beta,rr;
00988 double *d,*r,*p,*h;
00989
00990 d = new double [n];
00991 r = new double [n];
00992 p = new double [n];
00993 h = new double [n];
00994
00995
00996 if (iv==0){
00997 for (i=0;i<n;i++)
00998 x[i]=0.0;
00999 }
01000 mxv_scr (x,p);
01001 nory=0.0;
01002 for (i=0;i<n;i++){
01003 nory+=y[i]*y[i];
01004 r[i]=y[i]-p[i];
01005 }
01006
01007 if (nory<zero){
01008 print_err("right hand side vector has zero norm", __FILE__, __LINE__, __func__);
01009 ares=nory; ani=0;
01010 }
01011
01012 switch (tprec){
01013 case 1:{ prec_diag_scr (h,r); break; }
01014 case 5:{ prec_ssor_scr (h,r,par); break; }
01015 case 10:{ prec_ildl_scr (h,r); break; }
01016 case sparseindec:{
01017 sdirect->Solve (h,r);
01018 break;
01019 }
01020 default:{
01021 print_err("wrong preconditioner type is required", __FILE__, __LINE__, __func__);
01022 }
01023 }
01024
01025 nom = ss (r,h,n);
01026
01027 for (i=0;i<n;i++){
01028 d[i]=h[i];
01029 }
01030
01031
01032 for (i=0;i<ni;i++){
01033
01034
01035 mxv_scr (d,p);
01036
01037 denom = ss (d,p,n);
01038 if (fabs(denom)<zero){
01039 print_err("denominator in alpha expression is equal to zero", __FILE__, __LINE__, __func__);
01040 ares=nom; ani=i;
01041 return;
01042 }
01043
01044
01045 alpha = nom/denom;
01046
01047
01048 rr=0.0;
01049 for (j=0;j<n;j++){
01050 x[j]+=alpha*d[j];
01051 r[j]-=alpha*p[j];
01052 rr+=r[j]*r[j];
01053 }
01054
01055
01056
01057 switch (tprec){
01058 case 1:{ prec_diag_scr (h,r); break; }
01059 case 5:{ prec_ssor_scr (h,r,par); break; }
01060 case 10:{ prec_ildl_scr (h,r); break; }
01061 case sparseindec:{
01062 sdirect->Solve (h,r);
01063 break;
01064 }
01065 default:{
01066 print_err("wrong preconditioner type is required", __FILE__, __LINE__, __func__);
01067 }
01068 }
01069
01070 denom=nom;
01071 nom = ss (r,h,n);
01072 beta = nom/denom;
01073
01074 fprintf (stdout,"\n iteration %ld norres/norrhs %e",i,rr/nory);
01075
01076 if (fabs(rr/nory)<err){
01077 ani=i; ares=nom/nory;
01078 break;
01079 }
01080
01081
01082 for (j=0;j<n;j++){
01083 d[j]=beta*d[j]+h[j];
01084 }
01085 }
01086
01087 delete [] h;
01088 delete [] p;
01089 delete [] r;
01090 delete [] d;
01091 }
01092
01093
01094
01095
01096
01097
01098 void symcomprow::printmat (FILE *out)
01099 {
01100 long i,j;
01101
01102
01103
01104
01105
01106
01107
01108
01109
01110
01111
01112 fprintf (out,"\n\n\n matrix in the symmetric compressed row format\n");
01113 fprintf (out,"\n number of rows/columns %ld",n);
01114 fprintf (out,"\n number of stored entries %ld",negm);
01115 for (i=0;i<n;i++){
01116 fprintf (out,"\n row n. %ld ",i);
01117 for (j=adr[i];j<adr[i+1];j++){
01118 fprintf (out," (%ld, %f)",ci[j],a[j]);
01119 }
01120 }
01121 fprintf (out,"\n\n");
01122
01123 }
01124
01125
01126
01127
01128
01129
01130
01131
01132 void symcomprow::printdiag (FILE *out)
01133 {
01134 long i;
01135
01136 fprintf (out,"\n\n diagonal entries of matrix stored in the symmetric compressed row storage scheme\n");
01137
01138 for (i=0;i<n;i++){
01139 fprintf (out,"%lf\n",a[adr[i+1]-1]);
01140 }
01141 }
01142
01143
01144
01145
01146
01147
01148
01149
01150 void symcomprow::copy_scr (symcomprow &scr)
01151 {
01152 long i;
01153
01154 scr.n=n;
01155 scr.negm=negm;
01156
01157 if (scr.adr!=NULL){
01158 delete [] scr.adr;
01159 }
01160 scr.adr=new long [n+1];
01161
01162 if (scr.ci!=NULL){
01163 delete [] scr.ci;
01164 }
01165 scr.ci=new long [negm];
01166
01167 if (scr.a!=NULL){
01168 delete [] scr.a;
01169 }
01170 scr.a=new double [negm];
01171
01172 for (i=0;i<=n;i++){
01173 scr.adr[i]=adr[i];
01174 }
01175 for (i=0;i<negm;i++){
01176 scr.a[i]=a[i];
01177 scr.ci[i]=ci[i];
01178 }
01179 }
01180
01181
01182
01183
01184
01185
01186 long symcomprow::give_negm ()
01187 {
01188 return negm;
01189 }
01190
01191
01192
01193
01194
01195
01196
01197
01198
01199
01200
01201 void symcomprow::select_submatrix(symcomprow *smscr,long nsdof,long *sdof)
01202 {
01203
01204 long i,j,k,m,*auxci;
01205 long address1,address2;
01206
01207 smscr->n = nsdof;
01208
01209
01210
01211
01212 smscr->adr = new long[nsdof+1];
01213 auxci = new long[nsdof];
01214 for(i = 0; i < nsdof; i++){
01215 sdof[i]--;
01216 smscr->adr[i] = -1;
01217 auxci[i] = i;
01218 }
01219
01220
01221 smscr->negm = 0;
01222
01223 for(i = 0; i < nsdof; i++){
01224 address1 = adr[sdof[i]];
01225 address2 = adr[sdof[i]+1];
01226 m = 0;
01227 for(j = address1; j < address2; j++){
01228
01229 for(k = 0; k < nsdof; k++){
01230 if(ci[j] == sdof[k]){
01231 smscr->negm++;
01232 m++;
01233 break;
01234 }
01235 }
01236 }
01237
01238 smscr->adr[i] = smscr->negm-m;
01239 }
01240 smscr->adr[nsdof]= smscr->negm;
01241
01242
01243
01244 smscr->ci = new long[smscr->negm];
01245 smscr->a = new double[smscr->negm];
01246
01247 m = 0;
01248 for(i = 0; i < nsdof; i++){
01249 address1 = adr[sdof[i]];
01250 address2 = adr[sdof[i]+1];
01251 for(j = address1; j < address2; j++){
01252 for(k = 0; k < nsdof; k++){
01253 if(ci[j] == sdof[k]){
01254
01255 smscr->a[m] = a[j];
01256
01257 smscr->ci[m] = auxci[k];
01258 m++;
01259 }
01260 }
01261
01262 }
01263 }
01264
01265 delete []auxci;
01266
01267 }