00001 #include "dskyline.h"
00002 #include <math.h>
00003 #include <limits.h>
00004 #include <string.h>
00005
00006 dskyline::dskyline (void)
00007 {
00008
00009 n=0;
00010
00011 negm=0;
00012
00013
00014 adr = NULL;
00015
00016 a = NULL;
00017
00018 decompid=0;
00019 }
00020
00021
00022
00023 dskyline::~dskyline (void)
00024 {
00025 delete [] adr;
00026 delete [] a;
00027 }
00028
00029 long dskyline::decomp ()
00030 {
00031 return decompid;
00032 }
00033
00034 void dskyline::changedecomp ()
00035 {
00036 if (decompid==0) decompid=1;
00037 else decompid=0;
00038 }
00039
00040 void dskyline::setfact ()
00041 {
00042 decompid=1;
00043 }
00044 void dskyline::setnotfact ()
00045 {
00046 decompid=0;
00047 }
00048
00049
00050
00051
00052
00053
00054
00055
00056 void dskyline::allocadr (long m)
00057 {
00058 n=m;
00059 adr = new long [n+1];
00060 memset (adr,0,(n+1)*sizeof(long));
00061 }
00062
00063 double* dskyline::status ()
00064 {
00065 return a;
00066 }
00067
00068
00069
00070
00071
00072
00073
00074
00075
00076 void dskyline::column_lengths_elem (long *cn,long ndofe)
00077 {
00078 long i,j,k,min;
00079
00080 min=LONG_MAX;
00081 for (i=0;i<ndofe;i++){
00082 k=cn[i]-1;
00083 if (k>-1){
00084 if (k<min) min=k;
00085 }
00086 }
00087 for (i=0;i<ndofe;i++){
00088 k=cn[i]-1;
00089 if (k>-1){
00090 j=k-min+1;
00091 if (j>adr[k]) adr[k]=j;
00092 }
00093 }
00094 }
00095
00096
00097
00098
00099
00100
00101
00102
00103
00104
00105
00106 void dskyline::column_lengths_mult (long *ncn1,long *ncn2,long *mcn,long nm)
00107 {
00108 long i,j,k,min;
00109
00110 for (i=0;i<nm;i++){
00111
00112 if (ncn1[i]<ncn2[i])
00113 min=ncn1[i];
00114 else
00115 min=ncn2[i];
00116
00117
00118 j=mcn[i]-min+1;
00119
00120
00121 k=mcn[i]-1;
00122
00123
00124 if (j>adr[k])
00125 adr[k]=j;
00126 }
00127 }
00128
00129
00130
00131
00132
00133
00134
00135
00136
00137 void dskyline::column_lengths(gtopology *top)
00138 {
00139 long i,ne,ndofe,nm,*cn,*ncn1,*ncn2,*mcn;
00140
00141
00142 ne=top->ne;
00143
00144
00145 for (i=0;i<ne;i++){
00146 ndofe=top->give_ndofe (i);
00147 cn = new long [ndofe];
00148 top->give_code_numbers (i,cn);
00149 column_lengths_elem (cn,ndofe);
00150 delete [] cn;
00151 }
00152
00153
00154
00155
00156
00157 for (i=0;i<top->nen;i++){
00158
00159 nm=top->endnodes[i].ndofn;
00160 ncn1 = new long [nm];
00161 ncn2 = new long [nm];
00162 mcn = new long [nm];
00163
00164 top->give_endnode_code_numbers (i,ncn1,ncn2,mcn);
00165
00166 column_lengths_mult (ncn1,ncn2,mcn,nm);
00167 delete [] mcn;
00168 delete [] ncn2;
00169 delete [] ncn1;
00170 }
00171
00172
00173 for (i=0;i<top->nged;i++){
00174
00175 nm=top->gedges[i].ndofnf;
00176 ncn1 = new long [nm];
00177 ncn2 = new long [nm];
00178 mcn = new long [nm];
00179
00180 top->give_edge_code_numbers (i,1,ncn1,ncn2,mcn);
00181
00182 column_lengths_mult (ncn1,ncn2,mcn,nm);
00183 delete [] mcn;
00184 delete [] ncn2;
00185 delete [] ncn1;
00186
00187
00188 nm=top->gedges[i].ndofnl;
00189 ncn1 = new long [nm];
00190 ncn2 = new long [nm];
00191 mcn = new long [nm];
00192
00193 top->give_edge_code_numbers (i,2,ncn1,ncn2,mcn);
00194
00195 column_lengths_mult (ncn1,ncn2,mcn,nm);
00196 delete [] mcn;
00197 delete [] ncn2;
00198 delete [] ncn1;
00199 }
00200 }
00201
00202
00203
00204
00205
00206
00207
00208 void dskyline::addresses ()
00209 {
00210 long i,j,k;
00211
00212 j=adr[0]; adr[0]=0;
00213 for (i=1;i<n+1;i++){
00214 k=adr[i];
00215 adr[i]=j;
00216 j+=k;
00217 }
00218
00219 }
00220
00221
00222
00223
00224
00225
00226
00227
00228
00229 void dskyline::neglobmat ()
00230 {
00231 negm = 2*adr[n];
00232 }
00233
00234
00235
00236
00237
00238
00239 void dskyline::allocglomat ()
00240 {
00241 delete [] a;
00242 a = new double [negm];
00243 memset (a,0,negm*sizeof(double));
00244 }
00245
00246
00247
00248
00249
00250
00251 void dskyline::nullsky ()
00252 {
00253 long i;
00254 for (i=0;i<negm;i++){
00255 a[i]=0.0;
00256 }
00257 }
00258
00259
00260
00261
00262
00263
00264
00265
00266
00267
00268
00269 void dskyline::localize (matrix &b,long *cn)
00270 {
00271 long i,j,ii,jj,m;
00272
00273 m=adr[n];
00274
00275 for (i=0;i<b.m;i++){
00276 ii=cn[i]-1;
00277 if (ii<0) continue;
00278 for (j=0;j<b.m;j++){
00279 jj=cn[j]-1;
00280 if (jj<0) continue;
00281 if (ii<jj)
00282 a[adr[jj]+jj-ii+m] += b[i][j];
00283 else
00284 a[adr[ii]+ii-jj] += b[i][j];
00285 }
00286 }
00287 }
00288
00289
00290
00291
00292
00293
00294
00295
00296
00297
00298
00299
00300 void dskyline::localized (double *b,long *cn,long k)
00301 {
00302 long i,j,ii,jj,m;
00303
00304 m=adr[n];
00305
00306 for (i=0;i<k;i++){
00307 ii=cn[i]-1;
00308 if (ii==-1) continue;
00309 for (j=0;j<k;j++){
00310 jj=cn[j]-1;
00311 if (jj==-1) continue;
00312 if (ii<jj)
00313 a[adr[jj]+jj-ii+m] += b[i*k+j];
00314 else
00315 a[adr[ii]+ii-jj] += b[i*k+j];
00316 }
00317 }
00318 }
00319
00320
00321
00322
00323
00324
00325
00326
00327
00328
00329
00330 void dskyline::glocalize (matrix &b,long *rcn,long *ccn)
00331 {
00332 long i,j,ii,jj,m;
00333
00334 m=adr[n];
00335
00336 for (i=0;i<b.m;i++){
00337 ii=rcn[i]-1;
00338 if (ii<0) continue;
00339 for (j=0;j<b.n;j++){
00340 jj=ccn[j]-1;
00341 if (jj<0) continue;
00342 if (ii<jj)
00343 a[adr[jj]+jj-ii+m] += b[i][j];
00344 else
00345 a[adr[ii]+ii-jj] += b[i][j];
00346
00347 }
00348 }
00349
00350 }
00351
00352
00353
00354
00355
00356
00357
00358
00359
00360
00361
00362 void dskyline::mult_localize (long nm,long *ncn1,long *ncn2,long *mcn)
00363 {
00364 long i,ai,cn,m,mm;
00365
00366 m=adr[n];
00367
00368
00369 for (i=0;i<nm;i++){
00370 cn=mcn[i]-1;
00371 if (cn>-1){
00372
00373 ai=adr[cn];
00374
00375 mm=ncn1[i]-1;
00376 if (mm>-1){
00377 a[ai+cn-mm]=-1.0;
00378 a[ai+cn-mm+m]=-1.0;
00379 }
00380 mm=ncn2[i]-1;
00381 if (mm>-1){
00382 a[ai+cn-mm]=1.0;
00383 a[ai+cn-mm+m]=1.0;
00384 }
00385 }
00386 }
00387 }
00388
00389
00390 void dskyline::initiate (gtopology *top,long ndof,long mespr)
00391 {
00392 if (status()==NULL){
00393 allocadr (ndof);
00394 column_lengths (top);
00395 addresses ();
00396 neglobmat ();
00397 allocglomat ();
00398 }
00399 else{
00400 nullsky ();
00401 }
00402
00403 if (mespr==1) fprintf (stdout,"\n number of dskyline entries negm %ld",negm);
00404 }
00405
00406
00407
00408
00409
00410
00411
00412
00413
00414
00415 void dskyline::mxv_dsky (double *x,double *y)
00416 {
00417 long i,j,ii,m;
00418 double s;
00419
00420 m=adr[n];
00421
00422 for (i=0;i<n;i++){
00423 ii=i-(adr[i+1]-adr[i])+1; s=0.0;
00424 for (j=adr[i+1]-1;j>=adr[i];j--){
00425 s+=a[j]*x[ii]; ii++;
00426 }
00427 y[i]=s;
00428 }
00429 for (i=0;i<n;i++){
00430 ii=i-(adr[i+1]-adr[i])+1; s=x[i];
00431 for (j=adr[i+1]+m-1;j>adr[i]+m;j--){
00432 y[ii]+=a[j]*s; ii++;
00433 }
00434 }
00435 }
00436
00437
00438
00439
00440
00441
00442
00443
00444
00445 void dskyline::mtxv_dsky (double *x,double *y)
00446 {
00447 long i,j,ii,m;
00448 double s;
00449
00450 m=adr[n];
00451
00452 for (i=0;i<n;i++){
00453 ii=i-(adr[i+1]-adr[i])+1; s=0.0;
00454 for (j=adr[i+1]-1+m;j>=adr[i]+m;j--){
00455 s+=a[j]*x[ii]; ii++;
00456 }
00457 y[i]=s;
00458 }
00459 for (i=0;i<n;i++){
00460 ii=i-(adr[i+1]-adr[i])+1; s=x[i];
00461 for (j=adr[i+1]-1;j>=adr[i];j--){
00462 y[ii]+=a[j]*s; ii++;
00463 }
00464 }
00465 }
00466
00467 void dskyline::addmat_dsky (double c,dskyline &dsky)
00468 {
00469 long i;
00470 for (i=0;i<negm;i++){
00471 a[i]+=c*dsky.a[i];
00472 }
00473 }
00474
00475 void dskyline::scalmat_dsky (double c)
00476 {
00477 long i;
00478 for (i=0;i<negm;i++){
00479 a[i]*=c;
00480 }
00481 }
00482
00483
00484
00485
00486
00487
00488 void dskyline::copy_dsky (dskyline &dsky)
00489 {
00490 long i;
00491
00492
00493 dsky.n=n;
00494
00495 dsky.negm=negm;
00496
00497 if (dsky.adr!=NULL){
00498 delete [] dsky.adr;
00499 }
00500 dsky.adr = new long [n+1];
00501 if (dsky.a!=NULL){
00502 delete [] dsky.a;
00503 }
00504 dsky.a = new double [negm];
00505
00506 for (i=0;i<=n;i++){
00507 dsky.adr[i]=adr[i];
00508 }
00509 for (i=0;i<negm;i++){
00510 dsky.a[i]=a[i];
00511 }
00512 }
00513
00514
00515
00516
00517
00518
00519
00520
00521
00522
00523
00524
00525
00526
00527
00528
00529
00530
00531
00532
00533
00534
00535
00536
00537
00538
00539
00540 void dskyline::lu_dsky (double *x,double *y,double zero,long tc)
00541 {
00542 long i,j,k,ii,jj,kk,lj,lk,uk,m;
00543 double s;
00544
00545 m = adr[n];
00546
00547
00548
00549
00550 if (tc==1 || tc==2){
00551 if (decompid == 0){
00552
00553 for (i=0;i<n;i++){
00554
00555 lj=i-adr[i+1]+adr[i]+1;
00556
00557 for (j=lj;j<i;j++){
00558 kk=j-adr[j+1]+adr[j]+1;
00559 if (kk<lj) ii=lj;
00560 else ii=kk;
00561
00562 s=0.0; uk=adr[j]+j-ii; lk=adr[j]; jj=adr[i]+i-ii+m;
00563 for (k=uk;k>lk;k--){
00564 s+=a[k]*a[jj]; jj--;
00565 }
00566 a[jj]=(a[jj]-s)/a[lk];
00567
00568 s=0.0; uk+=m; lk+=m; jj=adr[i]+i-ii;
00569 for (k=uk;k>lk;k--){
00570 s+=a[jj]*a[k]; jj--;
00571 }
00572 a[jj]-=s;
00573 }
00574
00575 s=0.0; uk=adr[i+1]-1; lk=adr[i]; jj=adr[i+1]-1+m;
00576 for (k=uk;k>lk;k--){
00577 s+=a[jj]*a[k]; jj--;
00578 }
00579 a[lk]-=s;
00580 if (fabs(a[lk])<zero){
00581 fprintf (stderr,"\n\n zero diagonal entry in the dskyline (%ld-th row)",i);
00582 }
00583 }
00584
00585 decompid=1;
00586 }
00587 }
00588
00589
00590
00591
00592 if (tc==1 || tc==3){
00593 if (decompid == 0){
00594 print_err("matrix is not factorized, back substitution cannot be performed", __FILE__, __LINE__, __func__);
00595 }
00596 else{
00597 for (i=0;i<n;i++){
00598 ii=adr[i];
00599 s=0.0; lj=i+adr[i]-adr[i+1]+1; k=adr[i+1]-1;
00600 for (j=lj;j<i;j++){
00601 s+=a[k]*y[j]; k--;
00602 }
00603 if (fabs(a[ii])<zero){
00604 fprintf (stderr,"\n\n zero diagonal entry in the dskyline (%ld-th row)",i);
00605 }
00606 y[i]=(y[i]-s)/a[ii];
00607 }
00608
00609 for (i=n-1;i>-1;i--){
00610 x[i]=y[i]; s=x[i];
00611 lj=i-adr[i+1]+adr[i]+1; k=adr[i+1]-1+m;
00612 for (j=lj;j<i;j++){
00613 y[j]-=a[k]*s; k--;
00614 }
00615 }
00616 }
00617 }
00618 }
00619
00620
00621
00622
00623
00624
00625
00626
00627
00628
00629
00630
00631
00632
00633
00634 void dskyline::lukon_dsky (double *b,double *c,double *x,double *y,double zero,long nr,long tc)
00635 {
00636 long i,j,k,l,ii,jj,kk,lj,lk,uk,uj,m;
00637 double s;
00638
00639 m = adr[n];
00640
00641
00642
00643
00644 if (tc==1){
00645 if (decompid==0){
00646
00647 for (i=0;i<n-nr;i++){
00648
00649 lj=i-adr[i+1]+adr[i]+1;
00650
00651 for (j=lj;j<i;j++){
00652 kk=j-adr[j+1]+adr[j]+1;
00653 if (kk<lj) ii=lj;
00654 else ii=kk;
00655
00656 s=0.0; uk=adr[j]+j-ii; lk=adr[j]; jj=adr[i]+i-ii+m;
00657 for (k=uk;k>lk;k--){
00658 s+=a[k]*a[jj]; jj--;
00659 }
00660 a[jj]=(a[jj]-s)/a[lk];
00661
00662 s=0.0; uk+=m; lk+=m; jj=adr[i]+i-ii;
00663 for (k=uk;k>lk;k--){
00664 s+=a[jj]*a[k]; jj--;
00665 }
00666 a[jj]-=s;
00667 }
00668
00669 s=0.0; uk=adr[i+1]-1; lk=adr[i]; jj=adr[i+1]-1+m;
00670 for (k=uk;k>lk;k--){
00671 s+=a[jj]*a[k]; jj--;
00672 }
00673 a[lk]-=s;
00674 if (fabs(a[lk])<zero){
00675 fprintf (stderr,"\n\n zero diagonal entry (%e < %e) in the dskyline (%ld-th row), line %d",a[lk],zero,i,__LINE__);
00676 }
00677
00678
00679 }
00680
00681
00682 for (i=n-nr;i<n;i++){
00683 lj=i-adr[i+1]+adr[i]+1;
00684
00685 for (j=lj;j<n-nr;j++){
00686 kk=j-adr[j+1]+adr[j]+1;
00687 if (kk<lj) ii=lj;
00688 else ii=kk;
00689 s=0.0; uk=adr[j]+j-ii; lk=adr[j]; jj=adr[i]+i-ii+m;
00690 for (k=uk;k>lk;k--){
00691 s+=a[jj]*a[k]; jj--;
00692 }
00693 a[jj]=(a[jj]-s)/a[lk];
00694
00695 s=0.0; uk+=m; lk+=m; jj=adr[i]+i-ii;
00696 for (k=uk;k>lk;k--){
00697 s+=a[jj]*a[k]; jj--;
00698 }
00699 a[jj]-=s;
00700 }
00701
00702 for (j=n-nr;j<i;j++){
00703 kk=j-adr[j+1]+adr[j]+1;
00704 if (kk<lj) ii=lj;
00705 else ii=kk;
00706 s=0.0; uk=adr[j]+j-ii; lk=adr[j]+j-(n-nr); jj=adr[i]+i-ii+m;
00707 for (k=uk;k>lk;k--){
00708 s+=a[jj]*a[k]; jj--;
00709 }
00710 a[adr[i]+i-j+m]-=s;
00711
00712 s=0.0; uk=adr[i]+i-ii; lk=adr[i]+i-(n-nr); jj=adr[j]+j-ii+m;
00713 for (k=uk;k>lk;k--){
00714 s+=a[jj]*a[k]; jj--;
00715 }
00716 a[adr[i]+i-j]-=s;
00717
00718 }
00719
00720 s=0.0; uk=adr[i+1]-1; lk=adr[i]+i-(n-nr); jj=uk+m;
00721 for (k=uk;k>lk;k--){
00722 s+=a[k]*a[jj]; jj--;
00723 }
00724 a[adr[i]]-=s;
00725
00726
00727 }
00728
00729 decompid=1;
00730 }
00731
00732
00733
00734
00735
00736
00737
00738
00739
00740
00741
00742
00743
00744 long ack,ack1,acrk,ri,ci;
00745 for (i=n-nr;i<n;i++){
00746 ack=adr[i]; ack1=adr[i+1]-1;
00747 acrk=i-ack1+ack;
00748 if (acrk<n-nr){
00749 ri=i-n+nr; ci=i-n+nr;
00750 j=ack;
00751 b[ri*nr+ci]=a[j];
00752 ri--;
00753 for (j=ack+1;j<=i-n+nr+ack;j++){
00754 b[ri*nr+ci]=a[j+m];
00755 b[ci*nr+ri]=a[j];
00756
00757
00758 ri--;
00759 }
00760 }else{
00761 ri=i-n+nr; ci=i-n+nr;
00762 j=ack;
00763 b[ri*nr+ci]=a[j];
00764 ri--;
00765 for (j=ack+1;j<=ack1;j++){
00766 b[ri*nr+ci]=a[j+m];
00767 b[ci*nr+ri]=a[j];
00768
00769
00770 ri--;
00771 }
00772 }
00773 }
00774
00775
00776 for (i=0;i<n-nr;i++){
00777 ii=adr[i];
00778 s=0.0; lj=i+adr[i]-adr[i+1]+1; k=adr[i+1]-1;
00779 for (j=lj;j<i;j++){
00780 s+=a[k]*y[j]; k--;
00781 }
00782 if (fabs(a[ii])<zero){
00783 print_err("zero diagonal entry (%e < %e) in the dskyline (%ld-th row)",__FILE__, __LINE__, __func__,a[ii],zero,i);
00784 }
00785 y[i]=(y[i]-s)/a[ii];
00786 }
00787
00788 for (i=n-nr;i<n;i++){
00789 s=0.0; lj=i+adr[i]-adr[i+1]+1; k=adr[i+1]-1;
00790 for (j=lj;j<n-nr;j++){
00791 s+=a[k]*y[j]; k--;
00792 }
00793 y[i]-=s;
00794 }
00795
00796
00797 l=0;
00798 for (k=n-nr;k<n;k++){
00799 c[l]=y[k];
00800 l++;
00801 }
00802 }
00803
00804
00805
00806
00807 if (tc==2){
00808 if (decompid == 0){
00809 print_err("matrix is not factorized, back substitution cannot be performed", __FILE__, __LINE__, __func__);
00810 }
00811 else{
00812
00813
00814 j=0;
00815 for (k=n-nr;k<n;k++){
00816 x[k]=c[j]; j++;
00817
00818 }
00819
00820 for (i=n-1;i>=n-nr;i--){
00821 s=x[i];
00822 lj=adr[i]+1+m; uj=adr[i+1]+m; k=i-1;
00823 for (j=lj;j<uj;j++){
00824 y[k]-=a[j]*s; k--;
00825 }
00826 }
00827
00828 for (i=n-nr-1;i>-1;i--){
00829 x[i]=y[i]; s=x[i];
00830 lj=i-adr[i+1]+adr[i]+1; k=adr[i+1]-1+m;
00831 for (j=lj;j<i;j++){
00832 y[j]-=a[k]*s; k--;
00833 }
00834 }
00835 }
00836 }
00837
00838 }
00839
00840
00841
00842
00843
00844
00845
00846
00847
00848
00849
00850
00851
00852
00853
00854
00855
00856
00857
00858 void dskyline::bicg (double *x,double *y,long ni,double err,long &ani,double &ares,double zero,long iv)
00859 {
00860 long i,j;
00861 double normy,normr,normrr,nom,denom,alpha,beta;
00862 double *r,*rr,*d,*dd,*p,*pp;
00863
00864 r = new double [n]; rr = new double [n];
00865 d = new double [n]; dd = new double [n];
00866 p = new double [n]; pp = new double [n];
00867
00868 normy = ss (y,y,n);
00869
00870 if (normy<zero){
00871 fprintf (stderr,"\n\n norm of right hand side in bi-conjugate gradient method is smaller than %e",zero);
00872 fprintf (stderr,"\n see file %s, line %d.\n",__FILE__,__LINE__);
00873 ares=normy; ani=0;
00874 return;
00875 }
00876
00877 if (iv==0){
00878 for (i=0;i<n;i++)
00879 x[i]=0.0;
00880 }
00881
00882 mxv_dsky (x,p);
00883 for (i=0;i<n;i++){
00884 r[i]=y[i]-p[i];
00885 d[i]=r[i];
00886 rr[i]=r[i];
00887 dd[i]=rr[i];
00888 }
00889
00890 nom = ss (r,rr,n);
00891
00892 for (i=0;i<ni;i++){
00893 mxv_dsky (d,p);
00894 mtxv_dsky (dd,pp);
00895
00896 denom = ss (dd,p,n);
00897 alpha = nom/denom;
00898
00899 for (j=0;j<n;j++){
00900 x[j]+=alpha*d[j];
00901 r[j]-=alpha*p[j];
00902 rr[j]-=alpha*pp[j];
00903 }
00904
00905 normr = ss (r,r,n);
00906 normrr = ss (rr,rr,n);
00907 fprintf (stdout,"\n iteration %ld %e %e",i,normr/normy,normrr/normy);
00908 if (normr/normy<err && normrr/normy<err) break;
00909
00910
00911 denom = nom;
00912 nom = ss (r,rr,n);
00913 beta = nom/denom;
00914
00915 for (j=0;j<n;j++){
00916 d[j]=r[j]+beta*d[j];
00917 dd[j]=rr[j]+beta*dd[j];
00918 }
00919 }
00920
00921 ani=i; ares=normr;
00922 delete [] pp; delete [] p;
00923 delete [] dd; delete [] d;
00924 delete [] rr; delete [] r;
00925 }
00926
00927 void dskyline::printmat (FILE *out)
00928 {
00929 long i,j,m;
00930
00931
00932
00933
00934
00935
00936
00937
00938
00939
00940 m=adr[n];
00941
00942 fprintf (out,"\n\n");
00943 for (i=0;i<n;i++){
00944 fprintf (out,"\n");
00945 for (j=0;j<=i;j++){
00946 if (adr[i]+i-j>=adr[i+1]) fprintf (out," % 12.9le",0.0);
00947 else fprintf (out," % 12.9le",a[adr[i]+i-j]);
00948 }
00949 for (j=i+1;j<n;j++){
00950 if (adr[j]+j-i>=adr[j+1]) fprintf (out," % 12.9le",0.0);
00951 else fprintf (out," % 12.9le",a[adr[j]+j-i+m]);
00952 }
00953 }
00954 }
00955
00956 void dskyline::printdiag (FILE *out)
00957 {
00958 long i;
00959
00960 fprintf (out,"\n\n");
00961 for (i=0;i<n;i++){
00962 fprintf (out,"%5ld %20.10f\n",i,a[adr[i]]);
00963 }
00964 }
00965
00966
00967
00968
00969
00970
00971
00972 long dskyline::give_negm ()
00973 {
00974 return negm;
00975 }
00976
00977
00978
00979
00980
00981
00982
00983
00984
00985
00986
00987
00988
00989
00990 void dskyline::diag_check (double thr,double *rhs)
00991 {
00992 long i,j,nch;
00993
00994 nch=0;
00995 for (i=0;i<n;i++){
00996 j=adr[i];
00997 if (fabs(a[j])<thr){
00998 a[j]=1.0;
00999 rhs[i]=0.0;
01000 nch++;
01001 }
01002 }
01003 fprintf (stdout,"\n number of changes in the function dskyline::diag_check %ld\n",nch);
01004 }