00001 #include "locmatrix.h"
00002 #include <stdlib.h>
00003 #include <stdio.h>
00004
00005 locmatrix::locmatrix ()
00006 {
00007 nr=0; nc=0; nnc=0;
00008
00009 nncr=NULL; adr=NULL;
00010 lm=NULL; ci=NULL;
00011 }
00012
00013 locmatrix::~locmatrix ()
00014 {
00015 delete [] nncr;
00016 delete [] adr;
00017 delete [] lm;
00018 delete [] ci;
00019 }
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029 void locmatrix::initiate_var (long nrows,long ncolumns)
00030 {
00031 if (nrows!=nr){
00032 if (nncr!=NULL){
00033 delete [] nncr;
00034 nncr=NULL;
00035 }
00036 if (adr!=NULL){
00037 delete [] adr;
00038 adr=NULL;
00039 }
00040 }
00041
00042 nr=nrows;
00043 nc=ncolumns;
00044
00045 if (nncr==NULL){
00046 nncr = new long [nr];
00047 }
00048 }
00049
00050
00051
00052
00053
00054
00055 void locmatrix::addresses ()
00056 {
00057 long i;
00058
00059 if (adr==NULL){
00060 adr = new long [nr+1];
00061 }
00062
00063 adr[0]=0;
00064 for (i=1;i<nr+1;i++){
00065 adr[i]=adr[i-1]+nncr[i-1];
00066 }
00067
00068
00069 nnc=adr[nr];
00070 }
00071
00072
00073
00074
00075
00076
00077 void locmatrix::allocate_ci ()
00078 {
00079 if (ci!=NULL){
00080 delete [] ci;
00081 }
00082 ci = new long [nnc];
00083 }
00084
00085
00086
00087
00088
00089
00090 void locmatrix::allocate_lm ()
00091 {
00092 if (lm!=NULL){
00093 delete [] lm;
00094 }
00095 lm = new double [nnc];
00096 }
00097
00098
00099
00100
00101
00102
00103
00104
00105
00106 void locmatrix::initiate_nncr (long *a)
00107 {
00108 long i;
00109 for (i=0;i<nr;i++){
00110 nncr[i]=a[i];
00111 }
00112 }
00113
00114
00115
00116
00117
00118
00119
00120
00121
00122 void locmatrix::initiate_ci (long *colind)
00123 {
00124 long i;
00125 for (i=0;i<adr[nr];i++){
00126 ci[i]=colind[i];
00127 }
00128 }
00129
00130
00131
00132
00133
00134
00135
00136
00137
00138
00139 void locmatrix::initiate_ci (long id,long *colind)
00140 {
00141 long i,j;
00142
00143 if (id<0 || id>nr){
00144 fprintf (stderr,"\n wrong number of row is used in localization matrix (file %s, line %d)\n",__FILE__,__LINE__);
00145 }
00146
00147 j=0;
00148 for (i=adr[id];i<adr[id+1];i++){
00149 ci[i]=colind[i];
00150 j++;
00151 }
00152 }
00153
00154
00155
00156
00157
00158
00159
00160
00161
00162 void locmatrix::initiate_lm (double *a)
00163 {
00164 long i;
00165 for (i=0;i<adr[nr];i++){
00166 lm[i]=a[i];
00167 }
00168 }
00169
00170
00171
00172
00173
00174
00175
00176
00177
00178 void locmatrix::initiate_lm (long id,double *a)
00179 {
00180 long i,j;
00181
00182 if (id<0 || id>nr){
00183 fprintf (stderr,"\n wrong number of row is used in localization matrix (file %s, line %d)\n",__FILE__,__LINE__);
00184 }
00185
00186 j=0;
00187 for (i=adr[id];i<adr[id+1];i++){
00188 lm[i]=a[j];
00189 j++;
00190 }
00191 }
00192
00193
00194
00195
00196
00197
00198
00199
00200
00201 void locmatrix::lmxv (double *a,double *b)
00202 {
00203 long i,j,k;
00204 double s;
00205
00206 for (i=0;i<nr;i++){
00207 s=0.0;
00208 for (j=adr[i];j<adr[i+1];j++){
00209 k=ci[j];
00210 s+=lm[j]*a[k];
00211 }
00212 b[i]=s;
00213 }
00214 }
00215
00216
00217
00218
00219
00220
00221
00222 void locmatrix::lmtxv (double *a,double *b)
00223 {
00224 long i,j,k;
00225 double s;
00226
00227 for (i=0;i<nc;i++){
00228 b[i]=0.0;
00229 }
00230
00231 for (i=0;i<nr;i++){
00232 s=a[i];
00233 for (j=adr[i];j<adr[i+1];j++){
00234 k=ci[j];
00235 b[k]+=lm[j]*s;
00236 }
00237 }
00238 }
00239
00240
00241
00242
00243
00244
00245
00246
00247
00248
00249
00250 void locmatrix::lm01xv (double *a,double *b)
00251 {
00252 long i,j,k;
00253 double s;
00254
00255 for (i=0;i<nr;i++){
00256 s=0.0;
00257 for (j=adr[i];j<adr[i+1];j++){
00258 k=ci[j];
00259 s+=a[k];
00260 }
00261 b[i]=s;
00262 }
00263 }
00264
00265
00266
00267
00268
00269
00270
00271
00272
00273
00274
00275 void locmatrix::lmt01xv (double *a,double *b)
00276 {
00277 long i,j,k;
00278 double s;
00279
00280 for (i=0;i<nc;i++){
00281 b[i]=0.0;
00282 }
00283
00284 for (i=0;i<nr;i++){
00285 s=a[i];
00286 for (j=adr[i];j<adr[i+1];j++){
00287 k=ci[j];
00288 b[k]+=s;
00289 }
00290 }
00291 }
00292
00293
00294
00295
00296
00297
00298
00299
00300
00301
00302
00303
00304
00305
00306 void locmatrix::lm01xm (gmatrix &a,gmatrix &b)
00307 {
00308 long i,j,k,l;
00309 double s;
00310
00311 for (i=0;i<nr;i++){
00312 for (j=0;j<a.n;j++){
00313 s=0.0;
00314 for (k=adr[i];k<adr[i+1];k++){
00315 l=ci[k];
00316 s+=a.give_entry (i,l);
00317 }
00318 }
00319 b.add_entry (s,i,j);
00320 }
00321 }
00322
00323
00324
00325
00326
00327
00328
00329
00330
00331
00332
00333
00334
00335
00336
00337
00338 void locmatrix::lmxmxlmt01 (gmatrix &a,gmatrix &b)
00339 {
00340 long i,j,ii,jj,ncontr;
00341 long *rind,*cind,*aux;
00342 double s;
00343 double *entr;
00344
00345
00346 ncontr=0;
00347 for (i=0;i<nr;i++){
00348 ii=ci[i];
00349 for (j=0;j<nr;j++){
00350 jj=ci[j];
00351 s=a.give_entry (ii,jj);
00352 if (fabs(s)>threshold)
00353 ncontr++;
00354 }
00355 }
00356
00357
00358 rind = new long [ncontr];
00359 cind = new long [ncontr];
00360 entr = new double [ncontr];
00361
00362 ncontr=0;
00363 for (i=0;i<nr;i++){
00364 ii=ci[i];
00365 for (j=0;j<nr;j++){
00366 jj=ci[j];
00367 s=a.give_entry (ii,jj);
00368 if (fabs(s)>threshold){
00369 rind[ncontr]=ii;
00370 cind[ncontr]=jj;
00371 entr[ncontr]=s;
00372 ncontr++;
00373 }
00374 }
00375 }
00376
00377
00378
00379 aux = new long [nr+1];
00380 for (i=0;i<nr+1;i++){
00381 aux[i]=-1;
00382 }
00383
00384 for (i=0;i<ncontr;i++){
00385 if (rind[i]>cind[i]){
00386
00387 continue;
00388 }
00389 if (cind[i]-rind[i]+1>aux[cind[i]])
00390 aux[cind[i]]=cind[i]-rind[i]+1;
00391 }
00392
00393 ii=aux[0];
00394 aux[0]=0;
00395 for (i=1;i<nr+1;i++){
00396 jj=aux[i];
00397 aux[i]=aux[i-1]+ii;
00398 ii=jj;
00399 }
00400
00401
00402 b.sky->allocadr (nr);
00403
00404 for (i=0;i<nr+1;i++){
00405 b.sky->adr[i]=aux[i];
00406 }
00407
00408 b.sky->neglobmat ();
00409
00410 b.sky->allocglomat ();
00411
00412
00413 for (i=0;i<ncontr;i++){
00414 if (rind[i]>cind[i]){
00415
00416 continue;
00417 }
00418 j=aux[cind[i]];
00419 b.sky->a[j+cind[i]-rind[i]]=entr[i];
00420 }
00421
00422
00423 delete [] aux;
00424 delete [] rind;
00425 delete [] cind;
00426 delete [] entr;
00427 }
00428
00429
00430
00431
00432
00433
00434
00435
00436
00437
00438
00439
00440
00441
00442
00443 void locmatrix::lmxmxlmt (gmatrix &a,gmatrix &b)
00444 {
00445 long i,j,k,l,n,cind;
00446 long *aux,*nncaux,*adraux,**ciaux;
00447 double s;
00448
00449 n=a.n;
00450
00451 aux = new long [n];
00452
00453
00454 nncaux = new long [nr];
00455
00456
00457 adraux = new long [nr+1];
00458
00459
00460 ciaux = new long* [nr];
00461
00462
00463 adraux[0]=0;
00464 for (i=0;i<nr;i++){
00465 nncaux[i]=0;
00466 for (j=0;j<n;j++){
00467 s=0.0; l=adr[i]; aux[j]=-1;
00468 for (k=0;k<nncr[i];k++){
00469 cind=ci[l];
00470 s+=lm[l]*a.give_entry (i,cind);
00471 l++;
00472 }
00473 if (fabs(s)>threshold){
00474 nncaux[i]++;
00475 aux[j]=j;
00476 }
00477 }
00478 adr[i+1]=adr[i]+nncaux[i];
00479 ciaux[i] = new long [nncaux[i]];
00480 k=0;
00481 for (j=0;j<n;j++){
00482 if (aux[j]>-1){
00483 ciaux[i][k]=aux[j];
00484 k++;
00485 }
00486 }
00487 }
00488
00489
00490 }