00001 #include "skyline.h"
00002 #include <math.h>
00003 #include <limits.h>
00004 #include <string.h>
00005 #include <stdlib.h>
00006
00007 #include "cr.h"
00008 #include "scr.h"
00009
00010 skyline::skyline (void)
00011 {
00012
00013 n=0;
00014
00015 negm=0;
00016
00017
00018 adr = NULL;
00019
00020 a = NULL;
00021
00022
00023 decompid=0;
00024 }
00025
00026
00027 skyline::~skyline (void)
00028 {
00029 delete [] adr;
00030 delete [] a;
00031 }
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041 long skyline::decomp ()
00042 {
00043 return decompid;
00044 }
00045
00046
00047
00048
00049
00050
00051 void skyline::changedecomp ()
00052 {
00053 if (decompid==0) decompid=1;
00054 else decompid=0;
00055 }
00056
00057
00058
00059
00060 void skyline::setfact ()
00061 {
00062 decompid=1;
00063 }
00064
00065
00066
00067
00068 void skyline::setnotfact ()
00069 {
00070 decompid=0;
00071 }
00072
00073
00074
00075
00076
00077
00078
00079
00080
00081 void skyline::allocadr (long m)
00082 {
00083 n=m;
00084 adr = new long [n+1];
00085 memset (adr,0,(n+1)*sizeof(long));
00086 }
00087
00088
00089
00090
00091
00092
00093 double* skyline::status ()
00094 {
00095 return a;
00096 }
00097
00098
00099
00100
00101
00102
00103
00104
00105
00106 void skyline::column_lengths_elem (long *cn,long ndofe)
00107 {
00108 long i,j,k,min;
00109
00110 min=LONG_MAX;
00111 for (i=0;i<ndofe;i++){
00112 k=cn[i]-1;
00113 if (k>-1){
00114 if (k<min) min=k;
00115 }
00116 }
00117 for (i=0;i<ndofe;i++){
00118 k=cn[i]-1;
00119 if (k>-1){
00120 j=k-min+1;
00121 if (j>adr[k]) adr[k]=j;
00122 }
00123 }
00124 }
00125
00126
00127
00128
00129
00130
00131
00132
00133
00134
00135
00136 void skyline::column_lengths_mult (long *ncn1,long *ncn2,long *mcn,long nm)
00137 {
00138 long i,j,k,min;
00139
00140 for (i=0;i<nm;i++){
00141
00142 if (ncn1[i]<ncn2[i])
00143 min=ncn1[i];
00144 else
00145 min=ncn2[i];
00146
00147
00148 j=mcn[i]-min+1;
00149
00150
00151 k=mcn[i]-1;
00152
00153
00154 if (j>adr[k])
00155 adr[k]=j;
00156 }
00157 }
00158
00159
00160
00161
00162
00163
00164
00165
00166
00167
00168 void skyline::column_lengths (gtopology *top)
00169 {
00170 long i,ne,ndofe,nm,*cn,*ncn1,*ncn2,*mcn;
00171
00172
00173 ne=top->ne;
00174
00175
00176 for (i=0;i<ne;i++){
00177 ndofe=top->give_gndofe (i);
00178 cn = new long [ndofe];
00179 top->give_gcode_numbers (i,cn);
00180 column_lengths_elem (cn,ndofe);
00181 delete [] cn;
00182 }
00183
00184
00185
00186
00187
00188
00189
00190 for (i=0;i<top->nen;i++){
00191
00192 nm=top->endnodes[i].ndofn;
00193 ncn1 = new long [nm];
00194 ncn2 = new long [nm];
00195 mcn = new long [nm];
00196
00197 top->give_endnode_code_numbers (i,ncn1,ncn2,mcn);
00198
00199 column_lengths_mult (ncn1,ncn2,mcn,nm);
00200 delete [] mcn;
00201 delete [] ncn2;
00202 delete [] ncn1;
00203 }
00204
00205
00206 for (i=0;i<top->nged;i++){
00207
00208 nm=top->gedges[i].ndofnf;
00209 ncn1 = new long [nm];
00210 ncn2 = new long [nm];
00211 mcn = new long [nm];
00212
00213 top->give_edge_code_numbers (i,1,ncn1,ncn2,mcn);
00214
00215 column_lengths_mult (ncn1,ncn2,mcn,nm);
00216 delete [] mcn;
00217 delete [] ncn2;
00218 delete [] ncn1;
00219
00220
00221 nm=top->gedges[i].ndofnl;
00222 ncn1 = new long [nm];
00223 ncn2 = new long [nm];
00224 mcn = new long [nm];
00225
00226 top->give_edge_code_numbers (i,2,ncn1,ncn2,mcn);
00227
00228 column_lengths_mult (ncn1,ncn2,mcn,nm);
00229 delete [] mcn;
00230 delete [] ncn2;
00231 delete [] ncn1;
00232 }
00233 }
00234
00235
00236
00237
00238
00239
00240
00241
00242
00243 void skyline::addresses ()
00244 {
00245 long i,j,k;
00246
00247 j=adr[0]; adr[0]=0;
00248 for (i=1;i<n+1;i++){
00249 k=adr[i];
00250 adr[i]=j;
00251 j+=k;
00252 }
00253
00254 }
00255
00256
00257
00258
00259
00260
00261
00262
00263
00264 void skyline::neglobmat ()
00265 {
00266 negm = adr[n];
00267 }
00268
00269
00270
00271
00272
00273
00274
00275
00276 void skyline::allocglomat ()
00277 {
00278 delete [] a;
00279 a = new double [negm];
00280 memset (a,0,negm*sizeof(double));
00281 }
00282
00283
00284
00285
00286
00287
00288
00289
00290
00291
00292
00293
00294
00295 void skyline::localize (matrix &b,long *cn)
00296 {
00297 long i,j,ii,jj,k;
00298
00299 for (i=0;i<b.m;i++){
00300 ii=cn[i]-1;
00301 if (ii<0) continue;
00302 for (j=0;j<b.m;j++){
00303 jj=cn[j]-1;
00304 if (jj<0) continue;
00305 if (ii>jj) continue;
00306 k=adr[jj]+jj-ii;
00307 a[k]+=b[i][j];
00308 }
00309 }
00310
00311
00312
00313
00314
00315
00316
00317
00318
00319
00320
00321
00322
00323
00324
00325
00326
00327
00328 }
00329
00330
00331
00332
00333
00334
00335
00336
00337
00338
00339
00340
00341
00342
00343
00344 void skyline::localized (double *b,long *cn,long m)
00345 {
00346 long i,j,ii,jj,k;
00347
00348 for (i=0;i<m;i++){
00349 ii=cn[i]-1;
00350 if (ii<0) continue;
00351 for (j=0;j<m;j++){
00352 jj=cn[j]-1;
00353 if (jj<0) continue;
00354 if (ii>jj) continue;
00355 k=adr[jj]+jj-ii;
00356 a[k]+=b[i*m+j];
00357 }
00358 }
00359
00360
00361
00362
00363
00364
00365
00366
00367
00368
00369
00370
00371
00372
00373
00374
00375
00376 }
00377
00378
00379
00380
00381
00382
00383
00384
00385
00386
00387
00388 void skyline::glocalize (matrix &b,long *rcn,long *ccn)
00389 {
00390 long i,j,ii,jj,k;
00391
00392 for (i=0;i<b.m;i++){
00393 ii=rcn[i]-1;
00394 if (ii<0) continue;
00395 for (j=0;j<b.n;j++){
00396 jj=ccn[j]-1;
00397 if (jj<0) continue;
00398 if (ii>jj) continue;
00399 k=adr[jj]+jj-ii;
00400 a[k]+=b[i][j];
00401 }
00402 }
00403
00404 }
00405
00406
00407
00408
00409
00410
00411
00412
00413
00414
00415
00416 void skyline::mult_localize (long nm,long *ncn1,long *ncn2,long *mcn)
00417 {
00418 long i,ai,cn,m;
00419
00420
00421 for (i=0;i<nm;i++){
00422 cn=mcn[i]-1;
00423 if (cn>-1){
00424
00425 ai=adr[cn];
00426
00427 m=ncn1[i]-1;
00428 if (m>-1)
00429 a[ai+cn-m]=-1.0;
00430 m=ncn2[i]-1;
00431 if (m>-1)
00432 a[ai+cn-m]=1.0;
00433 }
00434 }
00435 }
00436
00437
00438
00439
00440
00441
00442 void skyline::nullsky ()
00443 {
00444 long i;
00445 for (i=0;i<negm;i++){
00446 a[i]=0.0;
00447 }
00448 }
00449
00450
00451
00452
00453
00454
00455
00456
00457 void skyline::copy_sky (skyline &sky)
00458 {
00459 long i;
00460
00461
00462 sky.n=n;
00463
00464 sky.negm=negm;
00465
00466 if (sky.adr!=NULL){
00467 delete [] sky.adr;
00468 }
00469 sky.adr = new long [n+1];
00470 if (sky.a!=NULL){
00471 delete [] sky.a;
00472 }
00473 sky.a = new double [negm];
00474
00475 for (i=0;i<=n;i++){
00476 sky.adr[i]=adr[i];
00477 }
00478 for (i=0;i<negm;i++){
00479 sky.a[i]=a[i];
00480 }
00481 }
00482
00483
00484
00485
00486
00487
00488
00489
00490
00491
00492 void skyline::initiate (gtopology *top,long ndof,long mespr)
00493 {
00494 if (status()==NULL){
00495 allocadr (ndof);
00496 column_lengths (top);
00497 addresses ();
00498 neglobmat ();
00499 allocglomat ();
00500 }
00501 else{
00502 nullsky ();
00503 }
00504
00505 if (mespr==1){
00506 fprintf (stdout,"\n number of skyline entries negm %ld",negm);
00507 fflush (stdout);
00508 }
00509
00510 }
00511
00512
00513
00514
00515
00516
00517
00518
00519
00520
00521
00522
00523
00524
00525
00526 void skyline::column_lengths_cr (long *cradr,long *ci)
00527 {
00528 long i,j;
00529
00530 for (i=0;i<n;i++){
00531 j=cradr[i];
00532 adr[i]=i-ci[j]+1;
00533 }
00534 }
00535
00536
00537
00538
00539
00540
00541
00542
00543
00544
00545 void skyline::mat_entries (long *cradr,long *ci,double *cra)
00546 {
00547 long i,j,lj,uj,k;
00548
00549 for (i=0;i<n;i++){
00550 lj=cradr[i];
00551 uj=cradr[i+1]-1;
00552 for (j=uj;j>=lj;j--){
00553 if (ci[j]<=i){
00554 k=adr[i]+i-ci[j];
00555 a[k]=cra[j];
00556 }
00557 }
00558 }
00559 }
00560
00561
00562
00563
00564
00565
00566
00567
00568
00569 void skyline::assemble_from_cr (long crn,long *cradr,long *crci,double *cra)
00570 {
00571
00572 allocadr (crn);
00573
00574
00575 column_lengths_cr (cradr,crci);
00576
00577
00578 addresses ();
00579
00580
00581 neglobmat ();
00582
00583
00584 allocglomat ();
00585
00586
00587 mat_entries (cradr,crci,cra);
00588
00589 fprintf (stdout,"\n skyline for aggregate contains %8ld unknowns and stores %10ld entries",n,negm);
00590 }
00591
00592
00593
00594
00595
00596
00597
00598
00599
00600
00601
00602
00603 void skyline::column_lengths_scr (long *scradr,long *scrci,long *se)
00604 {
00605 long i,j,k,lk,uk;
00606 long *aux;
00607 long min,max;
00608
00609 max=0;
00610 for (i=0;i<n;i++){
00611 if (se[i]>max)
00612 max=se[i];
00613 }
00614 aux = new long [max];
00615 for (i=0;i<max;i++){
00616 aux[i]=-1;
00617 }
00618 for (i=0;i<n;i++){
00619 aux[se[i]-1]=i;
00620 }
00621
00622
00623 for (i=0;i<n;i++){
00624
00625 j=se[i]-1;
00626
00627 lk=scradr[j];
00628 uk=scradr[j+1];
00629 min=max;
00630 for (k=lk;k<uk;k++){
00631 if (aux[scrci[k]]!=-1){
00632 if (aux[scrci[k]]<min)
00633 min=aux[scrci[k]];
00634 }
00635 }
00636 adr[i]=aux[j]-min+1;
00637 }
00638
00639 delete [] aux;
00640
00641
00642
00643
00644
00645
00646
00647
00648
00649
00650
00651
00652
00653
00654
00655
00656
00657
00658
00659
00660
00661
00662
00663
00664
00665 }
00666
00667
00668
00669
00670
00671 void skyline::mat_entries_scr (long *scradr,long *scrci,double *scra,long *se)
00672 {
00673 long i,j,k,l,lj,uj,ii,jj;
00674 long *aux;
00675
00676
00677
00678
00679 aux = new long [n];
00680 aux[0]=0;
00681 for (i=0;i<n-1;i++){
00682 aux[i+1]=aux[i]+se[i+1]-se[i]-1;
00683 }
00684
00685 for (i=0;i<n;i++){
00686 k=se[i]-1;
00687 lj=scradr[k];
00688 uj=scradr[k+1]-1;
00689 for (j=uj;j>lj-1;j--){
00690 for (ii=0;ii<n;ii++){
00691 if (scrci[j]==se[ii]-1){
00692 jj=scrci[j]-aux[ii];
00693 break;
00694 }
00695 }
00696
00697 l=adr[i] + k - aux[i] - jj;
00698 a[l]=scra[j];
00699 }
00700 }
00701
00702 delete [] aux;
00703 }
00704
00705
00706
00707
00708
00709
00710
00711
00712
00713
00714
00715
00716 void skyline::assemble_from_scr (long *scradr,long *scrci,double *scra,long neq,long *se)
00717 {
00718
00719 allocadr (neq);
00720
00721
00722 column_lengths_scr (scradr,scrci,se);
00723
00724
00725 addresses ();
00726
00727
00728 neglobmat ();
00729
00730
00731 allocglomat ();
00732
00733
00734 mat_entries_scr (scradr,scrci,scra,se);
00735
00736 fprintf (stdout,"\n skyline for subdomain/aggregate contains %8ld unknowns and stores %10ld entries",n,negm);
00737 }
00738
00739
00740
00741
00742
00743
00744
00745
00746
00747
00748
00749
00750
00751
00752
00753
00754
00755
00756 void skyline::ldl_sky (double *x,double *y,double zero,long tc)
00757 {
00758 long i,j,k,ac,ac1,ac2,acs,ack,ack1,acrk,aci,aci1,acri,acj,acj1;
00759 double s,g;
00760
00761
00762 if (tc==1 || tc==2){
00763 if (decompid == 0){
00764
00765
00766
00767 for (k=1;k<n;k++){
00768 ack=adr[k]; ack1=adr[k+1];
00769 ac1=k+ack; acrk=ac1-ack1+1;
00770 acj=ack1-2;
00771 for (i=acrk+1;i<k;i++){
00772 aci=adr[i]; aci1=adr[i+1];
00773 ac2=i+aci; acri=ac2-aci1+1;
00774 if (acri<acrk) ac=acrk;
00775 else ac=acri;
00776 acj1=ac1-ac; acs=ac2-ac;
00777 s=0.0;
00778 for (j=acj1;j>acj;j--){
00779 s+=a[j]*a[acs];
00780 acs--;
00781 }
00782 a[acj]-=s; acj--;
00783 }
00784 s=0.0;
00785 for (i=ack1-1;i>ack;i--){
00786 ac1=adr[acrk]; acrk++;
00787 g=a[i];
00788 a[i]/=a[ac1];
00789 s+=a[i]*g;
00790 }
00791 a[ack]-=s;
00792 if (fabs(a[ack])<zero){
00793 printf ("\n\n zero diagonal entry in the skyline (%ld-th row)",k);
00794 }
00795 }
00796
00797 decompid=1;
00798 }
00799 }
00800
00801 if (tc==1 || tc==3){
00802 if (decompid == 0){
00803 print_err("matrix is not factorized, back substitution cannot be performed", __FILE__, __LINE__, __func__);
00804 }
00805 else{
00806
00807
00808
00809
00810 for (k=1;k<n;k++){
00811 ack=adr[k]+1; ack1=adr[k+1];
00812 ac1=k-1; s=0.0;
00813 for (i=ack;i<ack1;i++){
00814 s+=a[i]*y[ac1];
00815 ac1--;
00816 }
00817 y[k]-=s;
00818 }
00819
00820 for (k=0;k<n;k++){
00821 ac=adr[k];
00822 if (fabs(a[ac])<zero){
00823 printf ("\n\n zero diagonal entry in the skyline (%ld-th row) (D^-1)",k);
00824 }
00825 y[k]/=a[ac];
00826 }
00827
00828 for (k=n-1;k>-1;k--){
00829 ack=adr[k]+1; ack1=adr[k+1];
00830 x[k]=y[k]; g=x[k];
00831 ac=k-1;
00832 for (i=ack;i<ack1;i++){
00833 y[ac]-=a[i]*g;
00834 ac--;
00835 }
00836 }
00837 }
00838 }
00839 }
00840
00841
00842
00843
00844
00845
00846
00847
00848
00849
00850
00851
00852
00853
00854
00855
00856
00857
00858
00859
00860
00861
00862
00863
00864
00865
00866
00867
00868
00869
00870
00871
00872
00873
00874
00875
00876
00877
00878
00879
00880
00881
00882
00883
00884
00885
00886
00887
00888
00889
00890
00891
00892
00893
00894
00895
00896
00897
00898
00899
00900
00901
00902
00903
00904
00905
00906
00907
00908
00909
00910
00911
00912
00913
00914
00915
00916
00917
00918
00919
00920
00921
00922
00923
00924
00925
00926
00927
00928
00929
00930
00931
00932
00933
00934
00935
00936
00937
00938
00939
00940
00941
00942
00943
00944
00945
00946
00947
00948
00949 long skyline::auxmax(long i,long j)
00950 {
00951 if (i>j) return i;
00952 else return j;
00953 }
00954
00955
00956 long skyline::auxmin(long i,long j)
00957 {
00958 if (j>i) return i;
00959 else return j;
00960 }
00961
00962
00963
00964
00965
00966
00967
00968
00969
00970
00971
00972
00973
00974
00975
00976
00977
00978
00979
00980
00981
00982
00983
00984
00985
00986
00987
00988
00989
00990
00991
00992
00993
00994
00995
00996
00997
00998
00999
01000
01001
01002
01003
01004
01005
01006
01007
01008
01009
01010
01011
01012
01013
01014
01015 void skyline::ldl_sky3 (double *x,double *y,long tc)
01016
01017
01018
01019
01020
01021
01022
01023
01024
01025
01026
01027
01028
01029
01030
01031
01032
01033
01034
01035
01036
01037
01038
01039 {
01040
01041 double sa,sb,s,g;
01042
01043 long i,j,k,ac1,ack,ack1,acrk,acj;
01044
01045 long aci,aci1,ac2,acri,ac,acj1,acs,pom;
01046
01047 long acia,aci1a,ac2a,acria,aca,acj1a,acsa,poma;
01048
01049 long acib,aci1b,ac2b,acrib,acb,acj1b,acsb,pomb,pombz;
01050
01051 long pocet;
01052
01053
01054
01055
01056
01057 if (tc==1 || tc==2){
01058
01059
01060
01061
01062
01063
01064
01065 for (k=1;k<n;k++){
01066
01067
01068
01069 if ((k%200)==0)
01070 {
01071
01072
01073
01074 }
01075 ack=adr[k];
01076
01077 ack1=adr[k+1];
01078
01079 ac1=k+ack;
01080
01081 acrk=ac1-ack1+1;
01082
01083 acj=ack1-2;
01084
01085
01086
01087 for (i=acrk+1;i<(k-1);i+=2){
01088
01089
01090
01091 acia=adr[i];
01092
01093 aci1a=acib=adr[i+1];
01094
01095 aci1b=adr[i+2];
01096
01097 ac2a=i+acia;
01098
01099 ac2b=i+1+acib;
01100
01101 acria=ac2a-aci1a+1;
01102
01103 acrib=ac2b-aci1b+1;
01104
01105 aca=auxmax(acrk,acria);
01106
01107 acb=auxmax(acrk,acrib);
01108
01109 acj1a=ac1-aca;
01110
01111 acj1b=ac1-acb;
01112
01113 acsa=ac2a-aca;
01114
01115 acsb=ac2b-acb;
01116
01117
01118
01119
01120
01121
01122
01123
01124
01125
01126
01127
01128
01129
01130
01131 poma=acsa+acj+1-acj1a;
01132
01133 pomb=pombz=acsb+acj+1-acj1b;
01134
01135 pocet=acj1a-acj;
01136
01137 #ifdef LADIC
01138
01139 if (pocet<0) pocet=0;
01140
01141 zvys_op(2l*pocet+1);
01142
01143 pocty[pocet]++;
01144
01145 #endif
01146
01147
01148
01149
01150
01151
01152
01153
01154
01155
01156
01157
01158
01159
01160
01161
01162
01163 sa=0.0;
01164
01165 sb=0.0;
01166
01167
01168
01169 if (acj1a<=acj1b)
01170
01171 {
01172
01173 for (j=acj+1;j<=acj1a;j++){
01174
01175 g=a[j];
01176
01177 sa+=g*a[poma++];
01178
01179 sb+=g*a[pomb++];
01180
01181 }
01182
01183 for (j=acj1a+1;j<=acj1b;j++){
01184
01185 sb+=a[j]*a[pomb++];
01186
01187 }
01188
01189 }
01190
01191 else
01192
01193 {
01194
01195 for (j=acj+1;j<=acj1b;j++){
01196
01197 g=a[j];
01198
01199 sa+=g*a[poma++];
01200
01201 sb+=g*a[pomb++];
01202
01203 }
01204
01205 for (j=acj1b+1;j<=acj1a;j++){
01206
01207 sa+=a[j]*a[poma++];
01208
01209 }
01210
01211 }
01212
01213
01214
01215 a[acj]-=sa;
01216
01217 sb+=a[acj]*a[pombz-1];
01218
01219 a[acj-1]-=sb;
01220
01221 acj-=2;
01222
01223
01224
01225 }
01226
01227 if (i<k)
01228
01229 {
01230
01231 aci=adr[i];
01232
01233 aci1=adr[i+1];
01234
01235 ac2=i+aci;
01236
01237 acri=ac2-aci1+1;
01238
01239 ac=auxmax(acrk,acri);
01240
01241
01242
01243 acj1=ac1-ac;
01244
01245 acs=ac2-ac;
01246
01247 s=0.0;
01248
01249 pom=acs+acj+1-acj1;
01250
01251 pocet=acj1-acj;
01252
01253 #ifdef LADIC
01254
01255 if (pocet<0) pocet=0;
01256
01257 zvys_op(2l*pocet+1);
01258
01259 pocty[pocet]++;
01260
01261 #endif
01262
01263
01264
01265
01266
01267 for (j=acj+1;j<=acj1;j++){
01268
01269 s+=a[j]*a[pom++];
01270
01271 }
01272
01273
01274
01275
01276
01277 a[acj]-=s; acj--;
01278
01279
01280
01281 }
01282
01283
01284
01285 s=0.0;
01286
01287 for (i=ack1-1;i>ack;i--){
01288
01289
01290
01291 ac1=adr[acrk]; acrk++;
01292
01293 g=a[i];
01294
01295 a[i]/=a[ac1];
01296
01297 s+=a[i]*g;
01298
01299 }
01300
01301
01302
01303 #ifdef LADIC
01304
01305 pocet=ack1-1-ack;
01306
01307 if (pocet<0) pocet=0;
01308
01309 zvys_op(3l*pocet+1l);
01310
01311 #endif
01312
01313 a[ack]-=s;
01314
01315 }
01316
01317 }
01318
01319 if (tc==1 || tc==3){
01320
01321
01322
01323
01324
01325
01326
01327
01328
01329 for (k=1;k<n;k++){
01330
01331
01332
01333 ack=adr[k]+1; ack1=adr[k+1];
01334
01335 ac1=k-1; s=0.0;
01336
01337 for (i=ack;i<ack1;i++){
01338
01339 s+=a[i]*y[ac1];
01340
01341 ac1--;
01342
01343 }
01344
01345 #ifdef LADIC
01346
01347 pocet=ack1-ack;
01348
01349 if (pocet<0) pocet=0;
01350
01351 zvys_op(2l*pocet+1l);
01352
01353 #endif
01354
01355 y[k]-=s;
01356
01357 }
01358
01359
01360
01361 for (k=0;k<n;k++){
01362
01363 ac=adr[k]; y[k]/=a[ac];
01364
01365 }
01366
01367 #ifdef LADIC
01368
01369 zvys_op(n);
01370
01371 #endif
01372
01373
01374
01375 for (k=n-1;k>-1;k--){
01376
01377
01378
01379 ack=adr[k]+1; ack1=adr[k+1];
01380
01381 x[k]=y[k]; g=x[k];
01382
01383 ac=k-1;
01384
01385 for (i=ack;i<ack1;i++){
01386
01387 y[ac]-=a[i]*g;
01388
01389 ac--;
01390
01391 }
01392
01393 #ifdef LADIC
01394
01395 pocet=ack1-ack;
01396
01397 if (pocet<0) pocet=0;
01398
01399 zvys_op(2l*pocet);
01400
01401 #endif
01402
01403 }
01404
01405 }
01406
01407 }
01408
01409 void skyline::eliminuj_4i_rev(double *a,long n,long s)
01410 {
01411 double sum1,sum2,sum3,sum4,x;
01412 long i,j,k;
01413
01414 for(j=0;j<n;j++)
01415 {
01416 sum1=sum2=0.0;
01417 for(k=0;k<(j-1);k+=2)
01418 {
01419 sum1+=(a[j*s+k]*a[j*s+k]);
01420 sum2+=(a[j*s+k+1]*a[j*s+k+1]);
01421 }
01422 if (j&1) sum1+=(a[j*s+j-1]*a[j*s+j-1]);
01423
01424 a[j*s+j]=sqrt(a[j*s+j]-sum1-sum2);
01425
01426 if(j&1)
01427 {
01428 for(i=j+1;i<(n-3);i+=4)
01429 {
01430 sum1=sum2=sum3=sum4=0.0;
01431 for(k=0;k<j;k++)
01432 {
01433 x=a[j*s+k];
01434 sum1+=a[i*s+k]*x;
01435 sum2+=a[(i+1)*s+k]*x;
01436 sum3+=a[(i+2)*s+k]*x;
01437 sum4+=a[(i+3)*s+k]*x;
01438 }
01439 x=a[j*s+j];
01440 a[i*s+j]=(a[i*s+j]-sum1)/x;
01441 a[(i+1)*s+j]=(a[(i+1)*s+j]-sum2)/x;
01442 a[(i+2)*s+j]=(a[(i+2)*s+j]-sum3)/x;
01443 a[(i+3)*s+j]=(a[(i+3)*s+j]-sum4)/x;
01444 }
01445 for(;i<n;i++)
01446 {
01447 sum1=0.0;
01448 for(k=0;k<j;k++)
01449 {
01450 sum1+=a[i*s+k]*a[j*s+k];
01451 }
01452 a[i*s+j]=(a[i*s+j]-sum1)/a[j*s+j];
01453 }
01454 }
01455 else
01456 {
01457 for(i=n-1;i>(j+3);i-=4)
01458 {
01459 sum1=sum2=sum3=sum4=0.0;
01460 for(k=0;k<j;k++)
01461 {
01462 x=a[j*s+k];
01463 sum1+=a[i*s+k]*x;
01464 sum2+=a[(i-1)*s+k]*x;
01465 sum3+=a[(i-2)*s+k]*x;
01466 sum4+=a[(i-3)*s+k]*x;
01467 }
01468 x=a[j*s+j];
01469 a[i*s+j]=(a[i*s+j]-sum1)/x;
01470 a[(i-1)*s+j]=(a[(i-1)*s+j]-sum2)/x;
01471 a[(i-2)*s+j]=(a[(i-2)*s+j]-sum3)/x;
01472 a[(i-3)*s+j]=(a[(i-3)*s+j]-sum4)/x;
01473 }
01474 for(;i>j;i--)
01475 {
01476 sum1=0.0;
01477 for(k=0;k<j;k++)
01478 {
01479 sum1+=a[i*s+k]*a[j*s+k];
01480 }
01481 a[i*s+j]=(a[i*s+j]-sum1)/a[j*s+j];
01482 }
01483 }
01484 }
01485 }
01486
01487
01488 void skyline::napln_a(long i1,long i2,long band,double *pole,double *a,long *adr)
01489 {
01490 long i,j,a1,a2,kolik;
01491 double *x;
01492 for(i=i1;i<i2;i++)
01493 {
01494 a1=adr[i];
01495 a2=adr[i+1];
01496 kolik=min(a2-a1,i+1-i1);
01497 x=pole+a1;
01498 for(j=0;j<kolik;j++) {a[(i-i1)*band+(i-i1)-j]=*x; x++; }
01499 for(j=kolik;j<=(i-i1);j++) a[(i-i1)*band+(i-i1)-j]=0.0;
01500 }
01501
01502 }
01503
01504
01505 void skyline::uloz_a(long i1,long i2,long band,double *pole,double *a,long *adr)
01506 {
01507 long i,j,a1,a2,kolik;
01508 double *x;
01509 for(i=i1;i<i2;i++)
01510 {
01511 a1=adr[i];
01512 a2=adr[i+1];
01513 kolik=min(a2-a1,i+1-i1);
01514 x=pole+adr[i];
01515 for(j=0;j<kolik;j++) { *x=a[(i-i1)*band+(i-i1)-j]; x++; }
01516 }
01517
01518 }
01519
01520
01521 void skyline::napln_b(long i1,long i2,long band,double *pole,double *b,long *adr)
01522 {
01523 long i,j,a1,a2,kolik;
01524 double *x;
01525 for(i=i1;i<i2;i++)
01526 {
01527 a1=adr[i];
01528 a2=adr[i+1];
01529 kolik=i+1-i1;
01530 x=pole+adr[i]+kolik;
01531 for(j=0;j<((a2-a1)-kolik);j++) { b[(i-i1)*band+(band-1)-j]=*x; x++; }
01532 for(j=((a2-a1)-kolik);j<band;j++) b[(i-i1)*band+(band-1)-j]=0.0;
01533 }
01534
01535 }
01536
01537
01538 void skyline::uloz_b(long i1,long i2,long band,double *pole,double *b,long *adr)
01539 {
01540 long i,j,a1,a2,kolik;
01541 double *x;
01542 for(i=i1;i<i2;i++)
01543 {
01544 a1=adr[i];
01545 a2=adr[i+1];
01546
01547 kolik=i+1-i1;
01548 x=pole+adr[i]+kolik;
01549 for(j=0;j<((a2-a1)-kolik);j++) { *x=b[(i-i1)*band+(band-1)-j]; x++; }
01550 }
01551
01552 }
01553
01554 void skyline::faze1_block3(double *a,double *b,long n,long s,long n2,long s2)
01555 {
01556
01557
01558 long i,j,k,i2,j2,stepi,stepj;
01559 double x,sum1,sum2,sum3,sum4;
01560
01561 stepi=1000;
01562 stepj=min(120000l/(16l*max(n,n2)),64l);
01563 for(i2=0;i2<n2;i2+=stepi)
01564 {
01565 for(j2=0;j2<n;j2+=stepj)
01566 {
01567 for(i=i2;i<(min(i2+stepi,n2)-3);i+=4)
01568 {
01569 for(j=j2;j<min(j2+stepj,n);j++)
01570 {
01571 sum1=sum2=sum3=sum4=0.0;
01572 for(k=0;k<j;k++)
01573 {
01574 x=a[j*s+k];
01575 sum1+=b[i*s2+k]*x;
01576 sum2+=b[(i+1)*s2+k]*x;
01577 sum3+=b[(i+2)*s2+k]*x;
01578 sum4+=b[(i+3)*s2+k]*x;
01579 }
01580 x=a[j*s+j];
01581 b[i*s2+j]=(b[i*s2+j]-sum1)/x;
01582 b[(i+1)*s2+j]=(b[(i+1)*s2+j]-sum2)/x;
01583 b[(i+2)*s2+j]=(b[(i+2)*s2+j]-sum3)/x;
01584 b[(i+3)*s2+j]=(b[(i+3)*s2+j]-sum4)/x;
01585 }
01586 }
01587 for(;i<min(i2+stepi,n2);i++)
01588 {
01589 for(j=j2;j<min(j2+stepj,n);j++)
01590 {
01591 sum1=0.0;
01592 for(k=0;k<j;k++) sum1+=b[i*s2+k]*a[j*s+k];
01593 b[i*s2+j]=(b[i*s2+j]-sum1)/a[j*s+j];
01594 }
01595 }
01596 }
01597 }
01598 }
01599
01600 void skyline::faze2_block3(double *a,double *b,long n,long s,long n2,long s2)
01601 {
01602 double x,sum1,sum2,sum3,sum4;
01603 long i,j,k,i2,j2,step;
01604
01605
01606
01607 step=120000l/(16l*max(n,n2));
01608 for(i2=0;i2<n2;i2+=step)
01609 {
01610 for(j2=0;j2<=i2;j2+=step)
01611 {
01612 for(i=i2;i<min(n2,i2+step);i++)
01613 {
01614 for(j=j2;j<=(min(i,j2+step-1)-3);j+=4)
01615 {
01616 sum1=sum2=sum3=sum4=0.0;
01617 for(k=0;k<n;k++)
01618 {
01619 x=b[i*s2+k];
01620 sum1+=x*b[j*s2+k];
01621 sum2+=x*b[(j+1)*s2+k];
01622 sum3+=x*b[(j+2)*s2+k];
01623 sum4+=x*b[(j+3)*s2+k];
01624 }
01625 a[i*s+j]-=sum1;
01626 a[i*s+j+1]-=sum2;
01627 a[i*s+j+2]-=sum3;
01628 a[i*s+j+3]-=sum4;
01629 }
01630 for(;j<=min(i,j2+step-1);j++)
01631 {
01632 sum1=0.0;
01633 for(k=0;k<n;k++) sum1+=b[i*s2+k]*b[j*s2+k];
01634 a[i*s+j]-=sum1;
01635 }
01636 }
01637 }
01638 }
01639 }
01640
01641 void skyline::blokove2(double *pole,long *adr,long n, long band)
01642 {
01643
01644
01645 long i,i1,i2;
01646 double *a,*b;
01647
01648
01649 a=new double [band*band];
01650 b=new double [band*band];
01651
01652 i=0;
01653
01654 i1=i2=band;
01655 napln_a(i,i1,band,pole,a,adr);
01656
01657
01658
01659
01660
01661
01662
01663
01664
01665
01666
01667 eliminuj_4i_rev(a,band,band);
01668
01669
01670
01671
01672
01673
01674
01675
01676 while(i1<n)
01677 {
01678
01679 i2=min(i1+band,n);
01680
01681
01682
01683 napln_b(i1,i2,i1-i,pole,b,adr);
01684
01685
01686
01687
01688
01689
01690
01691
01692
01693 faze1_block3(a,b,i1-i,i1-i,i2-i1,i1-i);
01694
01695
01696
01697
01698
01699
01700
01701
01702
01703
01704
01705
01706 uloz_a(i,i1,i1-i,pole,a,adr);
01707
01708
01709 napln_a(i1,i2,i2-i1,pole,a,adr);
01710
01711
01712
01713
01714
01715
01716
01717
01718
01719 faze2_block3(a,b,i1-i,i2-i1,i2-i1,i1-i);
01720
01721
01722
01723
01724
01725
01726
01727
01728
01729
01730
01731 uloz_b(i1,i2,i1-i,pole,b,adr);
01732
01733 eliminuj_4i_rev(a,i2-i1,i2-i1);
01734
01735
01736
01737
01738
01739
01740
01741
01742
01743 i=i1;
01744 i1=i2;
01745 }
01746
01747 uloz_a(i,n,i1-i,pole,a,adr);
01748
01749 }
01750
01751
01752
01753
01754
01755
01756
01757
01758
01759
01760
01761
01762
01763
01764
01765
01766
01767
01768
01769
01770
01771
01772
01773
01774
01775
01776
01777
01778
01779
01780
01781
01782
01783
01784
01785
01786
01787
01788
01789
01790
01791
01792
01793
01794
01795
01796
01797
01798
01799
01800
01801
01802
01803
01804
01805
01806
01807
01808
01809
01810
01811
01812
01813
01814
01815
01816
01817
01818
01819
01820
01821
01822
01823
01824
01825
01826
01827
01828
01829
01830
01831
01832
01833
01834
01835
01836
01837
01838
01839
01840
01841
01842
01843
01844 void skyline::mxv_sky (double *b,double *c)
01845 {
01846 long i,j,acb,aci,aci1;
01847 double s,g;
01848
01849 for (i=0;i<n;i++){
01850 aci=adr[i]; aci1=adr[i+1]-1;
01851 g=b[i]; s=0.0; acb=i-aci1+aci;
01852 for (j=aci1;j>aci;j--){
01853 s+=a[j]*b[acb];
01854 c[acb]+=a[j]*g;
01855 acb++;
01856 }
01857 c[i]=s+a[aci]*g;
01858 }
01859 }
01860
01861
01862
01863
01864
01865
01866
01867
01868
01869
01870
01871
01872
01873
01874
01875
01876
01877 void skyline::utv (double *b,double *c)
01878 {
01879 long i,j,k,lj,uj;
01880 double s;
01881
01882 nullv (c,n);
01883
01884 for (i=n-1;i>-1;i--){
01885 s=b[i]; k=i-1;
01886 c[i]+=b[i];
01887 lj=adr[i]+1; uj=adr[i+1];
01888 for (j=lj;j<uj;j++){
01889 c[k]+=a[j]*s; k--;
01890 }
01891 }
01892 }
01893
01894
01895
01896
01897
01898
01899
01900
01901
01902
01903
01904
01905
01906
01907
01908
01909
01910 void skyline::ltv (double *b,double *c)
01911 {
01912 long i,j,k,lj,uj;
01913 double s;
01914
01915 nullv (c,n);
01916
01917 for (i=0;i<n;i++){
01918 lj=adr[i]; uj=adr[i+1]-1;
01919 k=i-(uj-lj);
01920 s=0.0;
01921 for (j=uj;j>lj;j--){
01922 s+=a[j]*b[k]; k++;
01923 }
01924 c[i]=s+b[i];
01925 }
01926 }
01927
01928 void skyline::ldlmxv_sky (double *b,double *c)
01929 {
01930 long i;
01931
01932 utv (b,c);
01933 for (i=0;i<n;i++){
01934 c[i]*=a[adr[i]];
01935 }
01936 ltv (c,b);
01937
01938 for (i=0;i<n;i++){
01939 c[i]=b[i];
01940 }
01941 }
01942
01943
01944
01945
01946
01947
01948
01949
01950
01951
01952
01953 void skyline::addmat_sky (double c,skyline &sky)
01954 {
01955 long i;
01956 for (i=0;i<negm;i++){
01957 a[i]+=c*sky.a[i];
01958 }
01959 }
01960
01961
01962
01963
01964
01965
01966
01967
01968 void skyline::scalmat_sky (double c)
01969 {
01970 long i;
01971 for (i=0;i<negm;i++){
01972 a[i]*=c;
01973 }
01974 }
01975
01976
01977
01978
01979
01980
01981
01982
01983
01984
01985
01986
01987
01988
01989
01990
01991
01992
01993
01994
01995
01996
01997 void skyline::ker (double *r,long &nse,long *se,
01998 long ense,double limit,long tc)
01999 {
02000 long i,j,k,ii,jj,kk,lj,uj,li,ui,lk,uk,mi,ise,ib,*adrb;
02001 double s,g,*b;
02002
02003 b = new double [ense*n];
02004 adrb = new long [ense+1];
02005
02006
02007
02008
02009 if (tc==1 || tc==3){
02010 if (decompid == 1){
02011 print_err("factorized matrix has to be factorized once again", __FILE__, __LINE__, __func__);
02012 }
02013 else{
02014
02015
02016 ise=0;
02017 ib=0;
02018 adrb[0]=0;
02019
02020 for (i=1;i<n;i++){
02021 lj=adr[i]; uj=adr[i+1]-2;
02022
02023 mi=i-(uj-lj)-1; j=mi+1;
02024 for (jj=uj;jj>lj;jj--){
02025 li=adr[j]; ui=adr[j+1]-1; k=j-(ui-li);
02026
02027 if (k<mi) { uk=uj+1; ii=li+j-mi; }
02028 else { uk=lj+i-k; ii=ui; }
02029
02030 s=0.0;
02031 for (kk=uk;kk>jj;kk--){
02032 s+=a[kk]*a[ii]; ii--;
02033 }
02034
02035 a[jj]-=s; j++;
02036 }
02037
02038
02039 s=0.0; j=mi;
02040 for (jj=uj+1;jj>lj;jj--){
02041 g=a[jj];
02042 a[jj]/=a[adr[j]];
02043 s+=a[jj]*g;
02044 j++;
02045 }
02046 a[lj]-=s;
02047
02048
02049 if (fabs(a[lj])<limit){
02050
02051 se[ise]=i; ise++;
02052
02053 if (ise>ense){
02054 print_err("larger number of dependent rows than the estimated number", __FILE__, __LINE__, __func__);
02055 }
02056
02057
02058
02059 for (jj=uj+1;jj>lj;jj--){
02060 b[ib]=a[jj]; ib++;
02061 a[jj]=0.0;
02062 }
02063 a[lj]=1.0;
02064 adrb[ise]=ib;
02065
02066 for (j=i+1;j<n;j++){
02067 if (j-(adr[j+1]-adr[j])<i) a[adr[j]+j-i]=0.0;
02068 }
02069
02070 }
02071
02072 }
02073
02074 nse=ise;
02075
02076 decompid=1;
02077 }
02078 }
02079
02080
02081
02082
02083
02084 if (tc==2 || tc==3){
02085 if (tc==2 && decompid == 0){
02086 print_err("matrix is not factorized, base vectors of kernel cannot be computed", __FILE__, __LINE__, __func__);
02087 }
02088 else{
02089
02090 ise=nse;
02091 for (i=0;i<ise;i++){
02092 uj=adr[se[i]+1]-1; lj=adr[se[i]];
02093 ib=adrb[i];
02094 for (jj=uj;jj>lj;jj--){
02095 a[jj]=b[ib]; ib++;
02096 }
02097 }
02098
02099
02100 nullv (r,ise*n);
02101 for (i=0;i<ise;i++){
02102 ib=i*n;
02103 for (j=n-1;j>-1;j--){
02104 r[ib+se[i]]=1.0; s=r[ib+j];
02105 uk=adr[j+1]-1; lk=adr[j]; k=j-(uk-lk);
02106 for (kk=uk;kk>lk;kk--){
02107 r[ib+k]-=a[kk]*s; k++;
02108 }
02109 }
02110 }
02111
02112 }
02113 }
02114
02115 delete [] adrb; delete [] b;
02116 }
02117
02118
02119
02120
02121
02122
02123
02124
02125
02126
02127
02128
02129
02130
02131
02132
02133
02134
02135
02136
02137
02138
02139
02140
02141
02142
02143 void skyline::ldlkon_sky (double *b,double *c,double *x,double *y,long m,long tc)
02144 {
02145 long i,j,k,l,ac,ac1,aca,acc,acs,ack,ack1,acrk,ackk,aci,aci1,acri,acii;
02146 long acj,acj1,acui,acb1,acb2;
02147 double s,g;
02148
02149 acui=n-m;
02150 if (tc==1){
02151 if (decompid == 0){
02152
02153 for (k=1;k<acui;k++){
02154 ack=adr[k]; ack1=adr[k+1];
02155 ackk=k+ack; acrk=ackk-ack1+1;
02156 acj=ack1-2;
02157
02158 for (i=acrk+1;i<k;i++){
02159 aci=adr[i]; aci1=adr[i+1];
02160 acii=i+aci; acri=acii-aci1+1;
02161 if (acri<acrk) ac=acrk;
02162 else ac=acri;
02163 acj1=ackk-ac; acs=acii-ac;
02164 s=0.0;
02165 for (j=acj1;j>acj;j--){
02166 s+=a[j]*a[acs];
02167 acs--;
02168 }
02169 a[acj]-=s; acj--;
02170 }
02171
02172 s=0.0;
02173 for (i=ack1-1;i>ack;i--){
02174 ackk=adr[acrk]; acrk++;
02175 g=a[i];
02176 a[i]/=a[ackk];
02177 s+=a[i]*g;
02178 }
02179 a[ack]-=s;
02180 }
02181
02182
02183 for (k=acui;k<n;k++){
02184 ack=adr[k]; ack1=adr[k+1]-1;
02185 ackk=k+ack; acrk=ackk-ack1;
02186 acj=ack1-1;
02187 if (acrk<acui){
02188 for (i=acrk+1;i<acui;i++){
02189
02190 aci=adr[i]; aci1=adr[i+1];
02191 acii=i+aci; acri=acii-aci1+1;
02192 if (acri<acrk) ac=acrk;
02193 else ac=acri;
02194 acj1=ackk-ac; acs=acii-ac;
02195 s=0.0;
02196 for (j=acj1;j>acj;j--){
02197 s+=a[j]*a[acs];
02198 acs--;
02199 }
02200 a[acj]-=s; acj--;
02201 }
02202 aca=ackk-acui;
02203 for (i=acui;i<k;i++){
02204
02205 aci=adr[i]; aci1=adr[i+1];
02206 acii=i+aci; acri=acii-aci1+1;
02207 if (acri<acrk) ac=acrk;
02208 else ac=acri;
02209 acj=ackk-ac; acj1=ackk-acui;
02210 acs=acii-ac; s=0.0;
02211 for (j=acj;j>acj1;j--){
02212 s+=a[j]*a[acs];
02213 acs--;
02214 }
02215 a[aca]-=s; aca--;
02216 }
02217 s=0.0;
02218 for (i=acrk;i<acui;i++){
02219 aci=adr[i];
02220 g=a[ack1];
02221 a[ack1]/=a[aci];
02222 s+=a[ack1]*g;
02223 ack1--;
02224 }
02225 a[ack]-=s;
02226 }
02227 }
02228
02229 decompid=1;
02230 }
02231
02232
02233 acc=0;
02234 for (k=acui;k<n;k++){
02235 ack=adr[k]; ack1=adr[k+1]-1;
02236 acrk=k+ack-ack1;
02237 if (acrk<acui) aci=acui;
02238 else aci=acrk;
02239 aci1=ack+k-aci;
02240 acb1=acc*m+acc; acb2=acb1;
02241 for (i=ack;i<=aci1;i++){
02242 b[acb1]=a[i]; b[acb2]=a[i];
02243 acb1-=m; acb2--;
02244 }
02245 acc++;
02246 }
02247
02248
02249
02250
02251 for (k=0;k<acui;k++){
02252 ack=adr[k]; ack1=adr[k+1]-1;
02253 ackk=k+ack; acrk=ackk-ack1;
02254 s=0.0;
02255 for (i=ack1;i>ack;i--){
02256 s+=a[i]*y[acrk]; acrk++;
02257 }
02258 y[k]-=s;
02259 }
02260 for (k=acui;k<n;k++){
02261 ack=adr[k]; ack1=adr[k+1]-1;
02262 ackk=k+ack; acrk=ackk-ack1;
02263 s=0.0; aci=ackk-acui;
02264 for (i=ack1;i>aci;i--){
02265 s+=a[i]*y[acrk]; acrk++;
02266 }
02267 y[k]-=s;
02268 }
02269
02270
02271 l=0;
02272 for (k=acui;k<n;k++){
02273 c[l]=y[k];
02274 l++;
02275 }
02276
02277 }
02278
02279
02280 if (tc==2){
02281 if (decompid == 0){
02282 print_err("matrix is not factorized, back substitution cannot be performed", __FILE__, __LINE__, __func__);
02283 }
02284 else{
02285
02286
02287 j=0;
02288 for (k=acui;k<n;k++){
02289 x[k]=c[j]; j++;
02290 }
02291
02292
02293
02294
02295 for (k=0;k<acui;k++){
02296 ack=adr[k];
02297 y[k]/=a[ack];
02298 }
02299
02300 for (k=n-1;k>=acui;k--){
02301 ack=adr[k]; ack1=adr[k+1]-1;
02302 ackk=k+ack; acrk=ackk-ack1;
02303 s=x[k];
02304 for (i=ack1;i>ack;i--){
02305 y[acrk]-=a[i]*s;
02306 acrk++;
02307 }
02308 }
02309
02310 for (k=acui-1;k>-1;k--){
02311 ack=adr[k]; ack1=adr[k+1]-1;
02312 ackk=k+ack; acrk=ackk-ack1;
02313 x[k]=y[k]; s=x[k];
02314 for (i=ack1;i>ack;i--){
02315 y[acrk]-=s*a[i]; acrk++;
02316 }
02317 }
02318 }
02319 }
02320
02321 if (tc==3){
02322 if (decompid == 0){
02323 print_err("matrix is not factorized, back substitution cannot be performed", __FILE__, __LINE__, __func__);
02324 }
02325 else{
02326
02327
02328
02329
02330 for (k=1;k<acui;k++){
02331 ack=adr[k]+1; ack1=adr[k+1];
02332 ac1=k-1; s=0.0;
02333 for (i=ack;i<ack1;i++){
02334 s+=a[i]*y[ac1];
02335 ac1--;
02336 }
02337 y[k]-=s;
02338 }
02339
02340
02341
02342 for (k=0;k<acui;k++){
02343 ac=adr[k];
02344
02345
02346
02347 y[k]/=a[ac];
02348 }
02349
02350 for (k=acui-1;k>-1;k--){
02351 ack=adr[k]+1; ack1=adr[k+1];
02352 x[k]=y[k]; g=x[k];
02353 ac=k-1;
02354 for (i=ack;i<ack1;i++){
02355 y[ac]-=a[i]*g;
02356 ac--;
02357 }
02358 }
02359 }
02360 }
02361
02362 if (tc==4){
02363 if (decompid == 0){
02364 print_err ("matrix is not factorized, forward reduction of the right hand side cannot be performed", __FILE__, __LINE__, __func__);
02365 }
02366
02367
02368
02369 for (k=0;k<acui;k++){
02370 ack=adr[k]; ack1=adr[k+1]-1;
02371 ackk=k+ack; acrk=ackk-ack1;
02372 s=0.0;
02373 for (i=ack1;i>ack;i--){
02374 s+=a[i]*y[acrk]; acrk++;
02375 }
02376 y[k]-=s;
02377 }
02378 for (k=acui;k<n;k++){
02379 ack=adr[k]; ack1=adr[k+1]-1;
02380 ackk=k+ack; acrk=ackk-ack1;
02381 s=0.0; aci=ackk-acui;
02382 for (i=ack1;i>aci;i--){
02383 s+=a[i]*y[acrk]; acrk++;
02384 }
02385 y[k]-=s;
02386 }
02387
02388
02389 l=0;
02390 for (k=acui;k<n;k++){
02391 c[l]=y[k];
02392 l++;
02393 }
02394
02395 }
02396
02397 }
02398
02399
02400
02401
02402
02403
02404
02405
02406
02407
02408
02409
02410
02411
02412
02413
02414
02415
02416
02417
02418
02419
02420
02421
02422
02423 double skyline::ldlkoncount_sky (double *b,double *c,double *x,double *y,long m,long tc)
02424 {
02425 long i,j,k,ac,aca,acs,ack,ack1,acrk,ackk,aci,aci1,acri,acii;
02426 long acj,acj1,acui;
02427 double s;
02428
02429 double no;
02430
02431 no=0.0;
02432
02433 acui=n-m;
02434 if (tc==1){
02435 if (decompid == 0){
02436
02437 for (k=1;k<acui;k++){
02438 ack=adr[k]; ack1=adr[k+1];
02439 ackk=k+ack; acrk=ackk-ack1+1;
02440 acj=ack1-2;
02441
02442 for (i=acrk+1;i<k;i++){
02443 aci=adr[i]; aci1=adr[i+1];
02444 acii=i+aci; acri=acii-aci1+1;
02445 if (acri<acrk) ac=acrk;
02446 else ac=acri;
02447 acj1=ackk-ac; acs=acii-ac;
02448 s=0.0;
02449 for (j=acj1;j>acj;j--){
02450 no+=1.0;
02451 acs--;
02452 }
02453 acj--;
02454 }
02455
02456 s=0.0;
02457 for (i=ack1-1;i>ack;i--){
02458 ackk=adr[acrk]; acrk++;
02459 no+=1.0;
02460 }
02461 }
02462
02463 fprintf (stdout,"\n number of operations after inner unknowns elimination %le",no);
02464
02465
02466 for (k=acui;k<n;k++){
02467 ack=adr[k]; ack1=adr[k+1]-1;
02468 ackk=k+ack; acrk=ackk-ack1;
02469 acj=ack1-1;
02470 if (acrk<acui){
02471 for (i=acrk+1;i<acui;i++){
02472
02473 aci=adr[i]; aci1=adr[i+1];
02474 acii=i+aci; acri=acii-aci1+1;
02475 if (acri<acrk) ac=acrk;
02476 else ac=acri;
02477 acj1=ackk-ac; acs=acii-ac;
02478 s=0.0;
02479 for (j=acj1;j>acj;j--){
02480 no+=1.0;
02481 acs--;
02482 }
02483 acj--;
02484 }
02485 aca=ackk-acui;
02486 for (i=acui;i<k;i++){
02487
02488 aci=adr[i]; aci1=adr[i+1];
02489 acii=i+aci; acri=acii-aci1+1;
02490 if (acri<acrk) ac=acrk;
02491 else ac=acri;
02492 acj=ackk-ac; acj1=ackk-acui;
02493 acs=acii-ac; s=0.0;
02494 for (j=acj;j>acj1;j--){
02495 no+=1.0;
02496 acs--;
02497 }
02498 aca--;
02499 }
02500 s=0.0;
02501 for (i=acrk;i<acui;i++){
02502 aci=adr[i];
02503 no+=1.0;
02504 ack1--;
02505 }
02506 }
02507 }
02508
02509 decompid=1;
02510 }
02511
02512 }
02513
02514 return no;
02515 }
02516
02517
02518
02519
02520
02521
02522
02523
02524
02525
02526
02527
02528
02529
02530
02531
02532 void skyline::ldl_feti_sky (double *x,double *y,long nse,long *se,double zero)
02533 {
02534 long i,j,k,ii,lj,uj;
02535 double s;
02536
02537
02538
02539 k=0;
02540 for (i=0;i<n;i++){
02541 lj=adr[i]; uj=adr[i+1]-1; ii=i-uj+lj; s=0.0;
02542 for (j=uj;j>lj;j--){
02543 s+=a[j]*y[ii]; ii++;
02544 }
02545
02546 if (se[k]==i){
02547 y[i]=0.0; k++;
02548 }
02549 else{
02550 y[i]-=s;
02551 }
02552 }
02553
02554
02555 for (i=0;i<n;i++){
02556 s=a[adr[i]];
02557 if (fabs(s)<zero) fprintf (stderr,"\n\n zero diagonal entry (row number %ld) in function skyline::ldl_feti_sky\n",i);
02558 y[i]/=s;
02559 }
02560
02561
02562 k=nse-1;
02563 for (i=n-1;i>=0;i--){
02564 lj=adr[i]; uj=adr[i+1];
02565 if (k>-1){
02566 if (se[k]==i){
02567
02568 y[i]=0.0; k--;
02569 }
02570 }
02571 x[i]=y[i];
02572 s=x[i]; ii=i-1;
02573 for (j=lj+1;j<uj;j++){
02574 y[ii]-=s*a[j]; ii--;
02575 }
02576 }
02577
02578 }
02579
02580
02581
02582
02583
02584
02585
02586
02587
02588
02589
02590
02591
02592 void skyline::ldl_a12block (double *block,long nrdof)
02593 {
02594 long i,j,ai,ai1,li,lj,ri,ci,minri;
02595
02596 nullv (block,nrdof*(n-nrdof));
02597
02598 li=n-nrdof;
02599
02600 for (i=li;i<n;i++){
02601 ai=adr[i]; ai1=adr[i+1];
02602 minri=ai1-ai-i+1;
02603 if (minri>=li)
02604 continue;
02605 lj=ai+i-li+1;
02606 ri=li-1;
02607 ci=i-li;
02608 for (j=lj;j<ai1;j++){
02609 block[ri*nrdof+ci]=a[j];
02610 ri--;
02611 }
02612 }
02613
02614 }
02615
02616
02617
02618
02619
02620
02621
02622
02623 void skyline::printmat (FILE *out)
02624 {
02625 long i;
02626
02627 FILE *mat;
02628 mat = fopen ("matice.txt","w");
02629 fprintf (mat,"%ld",n);
02630 fprintf (out,"\n\n");
02631 for (i=0;i<n+1;i++){
02632 fprintf (out,"\n %9ld",adr[i]);
02633 }
02634 fprintf (out,"\n\n");
02635 for (i=0;i<negm;i++){
02636 fprintf (out,"\n %20.10f",a[i]);
02637 }
02638 fclose (mat);
02639
02640
02641
02642
02643
02644
02645
02646
02647
02648
02649
02650
02651
02652
02653
02654
02655
02656
02657
02658
02659
02660
02661
02662
02663
02664
02665
02666
02667
02668
02669
02670
02671
02672
02673
02674
02675
02676
02677
02678
02679
02680
02681
02682
02683
02684
02685
02686
02687
02688
02689
02690
02691
02692
02693
02694
02695
02696
02697
02698
02699
02700
02701
02702
02703
02704
02705
02706 }
02707
02708
02709
02710
02711
02712
02713
02714
02715 void skyline::printdiag(FILE *out)
02716 {
02717 fprintf (out,"\n\n");
02718 for (long i=0;i<n;i++)
02719 fprintf(out, "%5ld %.10e\n", i,a[adr[i]]);
02720 }
02721
02722
02723
02724
02725
02726
02727
02728
02729 double skyline::give_entry (long ri,long ci)
02730 {
02731 long i;
02732 double e;
02733
02734 if (ci>ri){
02735 i=adr[ci];
02736 e = a[i+ci-ri];
02737 }
02738 else{
02739 i=adr[ri];
02740 e = a[i+ri-ci];
02741 }
02742
02743 return e;
02744 }
02745
02746
02747
02748
02749
02750
02751
02752
02753
02754 void skyline::add_entry (double e,long ri,long ci)
02755 {
02756 long i;
02757
02758 if (ci>ri){
02759 i=adr[ci]+ci-ri;
02760 a[i]+=e;
02761 }
02762 else{
02763 i=adr[ri]+ri-ci;
02764 a[i]+=e;
02765 }
02766 }
02767
02768
02769
02770
02771
02772
02773 long skyline::give_negm ()
02774 {
02775 return negm;
02776 }
02777
02778
02779
02780
02781
02782
02783 void skyline::diag_scale (double *d)
02784 {
02785 long i,j,k;
02786
02787 for (i=0;i<n;i++){
02788 d[i]=1.0/sqrt(a[adr[i]]);
02789 }
02790
02791 for (i=0;i<n;i++){
02792 k=i;
02793 for (j=adr[i];j<adr[i+1];j++){
02794 a[j]*=d[i]*d[k];
02795 k--;
02796 }
02797 }
02798 }
02799
02800
02801
02802
02803
02804 void skyline::select_submatrix(skyline *smsky,long nsdof,long *sdof,FILE *out)
02805 {
02806 long i,j,k,l,address1,address2,upper,lower;
02807
02808 smsky->n = nsdof;
02809 smsky->adr = new long[nsdof+1];
02810 fprintf(out,"kontrola v select_submatrix\n");
02811
02812 for(i = 0; i < nsdof; i++){
02813 fprintf(out,"%ld\n",sdof[i]);
02814 smsky->adr[i]=-1;
02815 }
02816
02817
02818
02819
02820
02821 fprintf(out,"negm = %ld n = %ld\n",negm,n);
02822 for(i = 0; i < negm; i++){
02823 fprintf(out,"%le\n",a[i]);
02824 }
02825
02826 for(i = 0; i < n+1; i++){
02827 fprintf(out,"%ld\n",adr[i]);
02828 }
02829
02830
02831
02832 smsky->negm = 0;
02833 for(i = 0; i < nsdof; i++){
02834 address1 = adr[sdof[i]];
02835 address2 = adr[sdof[i]+1];
02836 k = address2 - address1;
02837 upper = sdof[i]+1;
02838 lower = upper-k;
02839 l = 0;
02840
02841 for(j = 0; j < nsdof; j++){
02842 if(sdof[j] < upper && sdof[j] >= lower){
02843
02844 smsky->negm++;
02845 l++;
02846 }
02847 }
02848 smsky->adr[i]=smsky->negm-l;
02849 }
02850
02851 fprintf(out,"smsky->negm = %ld\n",smsky->negm);
02852 smsky->adr[nsdof]= smsky->negm;
02853 for(i = 0; i < nsdof+1; i++){
02854 fprintf(out,"%ld\n",smsky->adr[i]);
02855 }
02856
02857 smsky->a = new double[smsky->negm];
02858
02859 for(i = 0; i < nsdof; i++){
02860 address1 = adr[sdof[i]];
02861 address2 = adr[sdof[i]+1];
02862 k = address2 - address1;
02863 upper = sdof[i]+1;
02864 lower = upper-k;
02865
02866 for(j = 0; j < nsdof; j++){
02867 if(sdof[j] < upper && sdof[j] >= lower){
02868
02869
02870
02871
02872 smsky->a[smsky->adr[i]+i-j]=a[address1+sdof[i]-sdof[j]];
02873
02874 }
02875 }
02876 }
02877
02878
02879
02880
02881 fprintf(out,"smsky negm = %ld n = %ld\n",smsky->negm,smsky->n);
02882 for(i = 0; i < smsky->negm; i++){
02883 fprintf(out,"%le\n",smsky->a[i]);
02884 }
02885
02886 for(i = 0; i < smsky->n+1; i++){
02887 fprintf(out,"%ld\n",smsky->adr[i]);
02888 }
02889
02890
02891
02892 }
02893
02894
02895
02896
02897
02898
02899
02900
02901
02902
02903
02904
02905
02906
02907 void skyline::diag_check (double thr)
02908 {
02909 long i,j;
02910
02911 for (i=0;i<n;i++){
02912 j=adr[i];
02913 if (fabs(a[j])<thr)
02914 a[j]=1.0;
02915 }
02916 }