00001
00002
00003 #include "SparseGridMtxLDL.h"
00004
00005 DSS_NAMESPASE_BEGIN
00006
00007
00008
00009 SparseGridMtxLDL::SparseGridMtxLDL(SparseMatrixF& sm,BYTE block_size,Ordering* block_order,MathTracer* eMT,BOOL load_data)
00010 : SparseGridMtx(sm,block_size,block_order,eMT)
00011 {
00012 IConectMatrix* bskl = block_order->cm;
00013 this->AlocateMemoryByPattern(bskl);
00014 ComputeBlocks();
00015
00016 if (sm.a!=NULL && load_data!=FALSE)
00017 LoadMatrixNumbers(sm);
00018 }
00019
00020
00021
00022 SparseGridMtxLDL::SparseGridMtxLDL(SparseMatrixF& sm,BYTE block_size,Ordering* block_order,Ordering* node_order,MathTracer* eMT,BOOL load_data )
00023 : SparseGridMtx(sm,block_size,block_order,node_order,eMT)
00024 {
00025 IConectMatrix* bskl = block_order->cm;
00026 this->AlocateMemoryByPattern(bskl);
00027 ComputeBlocks();
00028
00029 if (sm.a!=NULL && load_data!=FALSE)
00030 LoadMatrixNumbers(sm);
00031 }
00032
00033 void SparseGridMtxLDL::LoadZeros()
00034 {
00035 memset(this->Columns_data,0,columns_data_length*sizeof(double));
00036
00037 for (long bj=0; bj<n_blocks; bj++)
00038 for (long j=0; j<block_size; j++)
00039 Columns_data[bj*block_storage + block_size*j + j] = 1.0;
00040
00041 long* blockP = this->block_order->perm->Items;
00042 long* nodeP = this->node_order?node_order->perm->Items:NULL;
00043
00044 long neq = n_blocks*block_size - noDummyDOFs;
00045 for (long i=0; i<neq; i++)
00046 {
00047 long ni = (nodeP==NULL)?i:nodeP[i];
00048 long si = ni % block_size;
00049 long bi = ni / block_size;
00050 long nbi = (blockP==NULL)?bi:blockP[bi];
00051 Columns_data[nbi*block_storage+si+si*block_size] = 0.0;
00052 }
00053 }
00054
00055 void SparseGridMtxLDL::LoadMatrixNumbers(SparseMatrixF& sm)
00056 {
00057 LoadZeros();
00058 long* blockP = this->block_order->perm->Items;
00059 long* nodeP = this->node_order?node_order->perm->Items:NULL;
00060
00061 long aux_bi_idx = 0;
00062 long aux_bj_idx = 0;
00063
00064 for (ULONG j=0; j<sm.neq; j++)
00065 {
00066 long nj = (nodeP==NULL)?j:nodeP[j];
00067
00068 long sj = nj % block_size;
00069 long bj = nj / block_size;
00070 long nbj = (blockP==NULL)?bj:blockP[bj];
00071
00072 SparseGridColumn* Column_nbj = Columns[nbj];
00073 SparseGridColumn* Column_nbi = NULL;
00074
00075 double* column_nbj_bi_data = NULL;
00076 double* column_nbi_bj_data = NULL;
00077
00078 long old_nbi = -1;
00079
00080 for (ULONG ad = sm.Adr(j); ad<sm.Adr(j+1); ad++)
00081 {
00082 long i = (long)sm.Ci(ad);
00083 long ni = (nodeP==NULL)?i:nodeP[i];
00084 double val = sm.a[ad];
00085
00086 long si = ni % block_size;
00087 long bi = ni / block_size;
00088 long nbi = (blockP==NULL)?bi:blockP[bi];
00089
00090 if (old_nbi!=nbi)
00091 {
00092 Column_nbi = Columns[nbi];
00093 aux_bj_idx = aux_bi_idx = -1;
00094 if (nbi<nbj)
00095 {
00096 aux_bi_idx = Column_nbj->FindExistingBlockIndex(nbi);
00097 column_nbj_bi_data = Columns_data+Column_nbj->column_start_idx+aux_bi_idx*block_storage;
00098 }
00099 else if (nbi>nbj)
00100 {
00101 aux_bj_idx = Column_nbi->FindExistingBlockIndex(nbj);
00102 column_nbi_bj_data = Columns_data+Column_nbi->column_start_idx+aux_bj_idx*block_storage;
00103 }
00104 old_nbi=nbi;
00105 }
00106
00107 if (aux_bi_idx<0 && aux_bj_idx<0)
00108 {
00109 if (si<sj)
00110 Columns_data[nbi*block_storage+si+sj*block_size] = val;
00111 else
00112 Columns_data[nbi*block_storage+sj+si*block_size] = val;
00113 }
00114 else if (aux_bi_idx>=0)
00115 column_nbj_bi_data[block_size*sj + si] = val;
00116 else
00117 column_nbi_bj_data[block_size*si + sj] = val;
00118
00119
00120
00121 nonzeros++;
00122 }
00123 }
00124
00125
00126
00127
00128
00129
00130
00131
00132
00133
00134
00135
00136
00137
00138
00139
00140
00141
00142
00143
00144
00145
00146
00147
00148 }
00149
00150 SparseGridMtxLDL::~SparseGridMtxLDL()
00151 {
00152 if (Columns_data) {delete [] Columns_data; Columns_data = NULL;}
00153 }
00154
00155
00156 void SparseGridMtxLDL::AlocateMemoryByPattern(IConectMatrix* bskl)
00157 {
00158 clock_t ts = clock();
00159 this->n_blocks = bskl->N();
00160 columns_data_length = n_blocks*block_storage;
00161 SparseGridColumn* column;
00162
00163 for (long bj=0; bj<n_blocks; bj++)
00164 {
00165 IntArrayList* indexes = bskl->DetachIndexesAboveDiagonalInColumn(bj);
00166 Columns[bj] = column = new SparseGridColumn(indexes);
00167 column->column_start_idx = columns_data_length;
00168 columns_data_length += column->Entries*block_storage;
00169 }
00170
00171 if (columns_data_length==0)
00172 columns_data_length = 1;
00173
00174 long memn = ((long)columns_data_length)*sizeof(double)/1024;
00175 char str[256];
00176
00177 sprintf(str," sparse matrix size : %ld kB..",memn);
00178 Write(str);
00179 this->Columns_data = new double[columns_data_length];
00180 if (Columns_data==NULL)
00181 {
00182 Writeln("Out of memory!");
00183 return;
00184 }
00185
00186 clock_t end = clock();
00187 double duration = (double)(end - ts) / CLOCKS_PER_SEC;
00188 sprintf(str,"%.3f s",duration);
00189 Writeln(str);
00190
00191 }
00192
00193 double& SparseGridMtxLDL::ElementAt(int i, int j)
00194 {
00195 long* blockP = this->block_order->perm->Items;
00196 long* nodeP = this->node_order?node_order->perm->Items:NULL;
00197
00198 long aux_bi_idx = 0;
00199 long aux_bj_idx = 0;
00200
00201 long nj = (nodeP==NULL)?j:nodeP[j];
00202 long sj = nj % block_size;
00203 long bj = nj / block_size;
00204 long nbj = (blockP==NULL)?bj:blockP[bj];
00205
00206 SparseGridColumn* Column_nbj = Columns[nbj];
00207
00208 double* column_nbj_bi_data = NULL;
00209 double* column_nbi_bj_data = NULL;
00210
00211 long ni = (nodeP==NULL)?i:nodeP[i];
00212 long si = ni % block_size;
00213 long bi = ni / block_size;
00214 long nbi = (blockP==NULL)?bi:blockP[bi];
00215
00216 SparseGridColumn* Column_nbi = Columns[nbi];
00217
00218 aux_bj_idx = aux_bi_idx = -1;
00219 if (nbi<nbj)
00220 {
00221 aux_bi_idx = Column_nbj->FindExistingBlockIndex(nbi);
00222 column_nbj_bi_data = Columns_data+Column_nbj->column_start_idx+aux_bi_idx*block_storage;
00223 }
00224 else if (nbi>nbj)
00225 {
00226 aux_bj_idx = Column_nbi->FindExistingBlockIndex(nbj);
00227 column_nbi_bj_data = Columns_data+Column_nbi->column_start_idx+aux_bj_idx*block_storage;
00228 }
00229
00230 if (aux_bi_idx<0 && aux_bj_idx<0)
00231 {
00232 if (si<sj)
00233 return Columns_data[nbi*block_storage+si+sj*block_size];
00234 else
00235 return Columns_data[nbi*block_storage+sj+si*block_size];
00236 }
00237 else if (aux_bi_idx>=0)
00238 return column_nbj_bi_data[block_size*sj + si];
00239 else
00240 return column_nbi_bj_data[block_size*si + sj];
00241 }
00242
00243
00244
00245
00246
00247
00248
00249
00250
00251
00252
00253
00254
00255
00256
00257
00258
00259
00260
00261
00262
00263
00264
00265
00266
00267
00268
00269
00270
00271
00272
00273
00274
00275
00276
00277
00278
00279
00280 void SparseGridMtxLDL::MultiplyByVector(const LargeVectorAttach& x, LargeVectorAttach& y)
00281 {
00282 y.Initialize();
00283
00284 double* py = y.DataPtr();
00285 double* px = x.DataPtr();
00286 for (long bj=0; bj<n_blocks; bj++)
00287 {
00288 SparseGridColumn& columnJ = *Columns[bj];
00289
00290 long cnt = columnJ.Entries;
00291 for (long idx = 0; idx<cnt; idx++)
00292 {
00293
00294 long bi = columnJ.IndexesUfa->Items[idx];
00295
00296 long block = columnJ.column_start_idx + idx*block_storage;
00297 BlockArith->AddMultBlockByVectorSym(Columns_data+block,px,py,block_size*bi,bj*block_size);
00298 }
00299
00300 BlockArith->MultDiagonalBlockByVector(Columns_data+bj*block_storage,px+bj*block_size, py+bj*block_size);
00301 }
00302 }
00303
00304 void SparseGridMtxLDL::Factorize()
00305 {
00306
00307
00308
00309
00310 double* Atmp = new double[block_storage];
00311 long bi,bj,min_bi_J = 0;
00312
00313
00314 long* p_blockJ_pattern = new long[n_blocks+1];
00315 memset(p_blockJ_pattern,0,(n_blocks+1)*sizeof(long));
00316
00317 long newdotcnount = (long)ceil((double)n_blocks / 24.0);
00318
00319 BlockArith->zero_pivots = 0;
00320
00321 no_multiplications = 0;
00322
00323 double* cd = this->Columns_data;
00324 double* atmp = Atmp;
00325 {
00326 double* dd = cd;
00327 double* idd = cd;
00328 long Djj = 0;
00329 eMT->act_block = 0;
00330 for (bj=0; bj<n_blocks;bj++)
00331 {
00332 SparseGridColumn &columnJ = *Columns[bj];
00333 long noJentries = columnJ.Entries;
00334 if (noJentries>0)
00335 {
00336 long* columnJentries = columnJ.IndexesUfa->Items;
00337 double* pAkj = cd+columnJ.column_start_idx;
00338 double* pAij = pAkj;
00339
00340
00341 min_bi_J=bj;
00342 for (long i=noJentries-1; i>=0; i--)
00343 p_blockJ_pattern[min_bi_J=columnJentries[i]] = ~(i*block_storage);
00344
00345
00346 for (long idx_J=1; idx_J<noJentries; idx_J++)
00347 {
00348 pAij += block_storage;
00349 bi = columnJentries[idx_J];
00350
00351 SparseGridColumn &columnI = *Columns[bi];
00352 long noIentries = columnI.Entries;
00353
00354 if (noIentries>0)
00355 {
00356 double* pAki = cd+columnI.column_start_idx + (noIentries-1)*block_storage;
00357 long* columnIentries = columnI.IndexesUfa->Items;
00358 for (long* columnIentry = columnIentries+noIentries-1; columnIentry>=columnIentries; pAki-=block_storage)
00359 {
00360 long idx_K = p_blockJ_pattern[*columnIentry--];
00361 if ( idx_K==0) continue;
00362 BlockArith->SubATBproduct(pAij,pAki,pAkj+ ~idx_K );
00363
00364 }
00365 }
00366 }
00367
00368
00369
00370 for (long idx=noJentries-1; idx>=0; idx--)
00371 {
00372 bi = columnJentries[idx];
00373
00374 p_blockJ_pattern[bi] = 0;
00375
00376
00377 long Aij = columnJ.column_start_idx + idx*block_storage;
00378
00379
00380 Array::Copy(this->Columns_data,Aij,Atmp,0,block_storage);
00381
00382
00383 BlockArith->SubstSolveBlock(idd+bi*block_storage,cd+Aij);
00384
00385
00386
00387
00388 BlockArith->SubATBproduct(dd+Djj,atmp,cd+Aij);
00389 }
00390 }
00391
00392
00393 BlockArith->FactorizeBlock(dd+Djj);
00394
00395 if ((bj % newdotcnount) == 0)
00396 Write(".");
00397
00398 Djj += block_storage;
00399 eMT->act_block += block_size;
00400 if(eMT->break_flag)
00401 break;
00402 }
00403 }
00404 delete [] p_blockJ_pattern;
00405 delete [] Atmp;
00406 ComputeBlocks();
00407 }
00408
00409 void SparseGridMtxLDL::FactorizeOMP()
00410 {
00411
00412
00413 double* Atmp = new double[block_storage];
00414 long bi,bj,min_bi_J = 0;
00415
00416
00417 long* p_blockJ_pattern = new long[n_blocks+1];
00418 memset(p_blockJ_pattern,0,(n_blocks+1)*sizeof(long));
00419
00420 long newdotcnount = (long)ceil((double)n_blocks / 24.0);
00421
00422 BlockArith->zero_pivots = 0;
00423
00424 no_multiplications = 0;
00425
00426 double* cd = this->Columns_data;
00427 double* atmp = Atmp;
00428 {
00429 double* dd = cd;
00430 double* idd = cd;
00431 long Djj = 0;
00432 eMT->act_block = 0;
00433 for (bj=0; bj<n_blocks;bj++)
00434 {
00435 SparseGridColumn &columnJ = *Columns[bj];
00436 long noJentries = columnJ.Entries;
00437 if (noJentries>0)
00438 {
00439 long* columnJentries = columnJ.IndexesUfa->Items;
00440 double* pAkj = cd+columnJ.column_start_idx;
00441 double* pAij = pAkj;
00442
00443
00444 min_bi_J=bj;
00445 for (long i=noJentries-1; i>=0; i--)
00446 p_blockJ_pattern[min_bi_J=columnJentries[i]] = ~(i*block_storage);
00447
00448
00449 int chunk = 10000;
00450
00451 #pragma omp parallel
00452 {
00453 #pragma omp for schedule(dynamic,chunk)
00454 for (long idx_J=1; idx_J<noJentries; idx_J++)
00455 {
00456 pAij += block_storage;
00457 bi = columnJentries[idx_J];
00458
00459 SparseGridColumn &columnI = *Columns[bi];
00460 long noIentries = columnI.Entries;
00461
00462 if (noIentries>0)
00463 {
00464 double* pAki = cd+columnI.column_start_idx + (noIentries-1)*block_storage;
00465 long* columnIentries = columnI.IndexesUfa->Items;
00466 for (long* columnIentry = columnIentries+noIentries-1; columnIentry>=columnIentries; pAki-=block_storage)
00467 {
00468 long idx_K = p_blockJ_pattern[*columnIentry--];
00469 if ( idx_K==0) continue;
00470 BlockArith->SubATBproduct(pAij,pAki,pAkj+ ~idx_K );
00471
00472 }
00473 }
00474 }
00475 }
00476
00477
00478
00479 for (long idx=noJentries-1; idx>=0; idx--)
00480 {
00481 bi = columnJentries[idx];
00482
00483 p_blockJ_pattern[bi] = 0;
00484
00485
00486 long Aij = columnJ.column_start_idx + idx*block_storage;
00487
00488
00489 Array::Copy(this->Columns_data,Aij,Atmp,0,block_storage);
00490
00491
00492 BlockArith->SubstSolveBlock(idd+bi*block_storage,cd+Aij);
00493
00494
00495
00496
00497 BlockArith->SubATBproduct(dd+Djj,atmp,cd+Aij);
00498 }
00499 }
00500
00501
00502 BlockArith->FactorizeBlock(dd+Djj);
00503
00504 if ((bj % newdotcnount) == 0)
00505 Write(".");
00506
00507 Djj += block_storage;
00508 eMT->act_block += block_size;
00509 if(eMT->break_flag)
00510 break;
00511 }
00512 }
00513 delete [] p_blockJ_pattern;
00514 delete [] Atmp;
00515 ComputeBlocks();
00516 }
00517
00518 void SparseGridMtxLDL::Factorize_Incomplete()
00519 {
00520 double* Atmp = new double[block_storage];
00521 long bi,bj,min_bi_J = 0;
00522
00523
00524 long* p_blockJ_pattern = new long[n_blocks+1];
00525 memset(p_blockJ_pattern,0,(n_blocks+1)*sizeof(long));
00526
00527 long newdotcnount = (long)ceil((double)n_blocks / 24.0);
00528
00529 BlockArith->zero_pivots = 0;
00530
00531 double* cd = this->Columns_data;
00532 double* atmp = Atmp;
00533 {
00534 double* dd = cd;
00535 double* idd = cd;
00536 long Djj = 0;
00537 eMT->act_block = 0;
00538 for (bj=0; bj<n_blocks;bj++)
00539 {
00540 SparseGridColumn &columnJ = *Columns[bj];
00541 long noJentries = columnJ.Entries;
00542 if (noJentries>0)
00543 {
00544 long* columnJentries = columnJ.IndexesUfa->Items;
00545 double* pAkj = cd+columnJ.column_start_idx;
00546 double* pAij = pAkj;
00547
00548
00549 min_bi_J=bj;
00550 for (long i=noJentries-1; i>=0; i--)
00551 p_blockJ_pattern[min_bi_J=columnJentries[i]] = ~(i*block_storage);
00552
00553
00554 for (long idx_J=1; idx_J<noJentries; idx_J++)
00555 {
00556 pAij += block_storage;
00557 bi = columnJentries[idx_J];
00558
00559
00560 SparseGridColumn &columnI = *Columns[bi];
00561 long noIentries = columnI.Entries;
00562
00563 if (noIentries>0)
00564 {
00565 double* pAki = cd+columnI.column_start_idx + (noIentries-1)*block_storage;
00566 long* columnIentries = columnI.IndexesUfa->Items;
00567 for (long* columnIentry = columnIentries+noIentries-1; columnIentry>=columnIentries; pAki-=block_storage)
00568 {
00569 long idx_K = p_blockJ_pattern[*columnIentry--];
00570 if (idx_K==0) continue;
00571 BlockArith->SubATBproduct(pAij,pAki,pAkj+ ~idx_K );
00572 }
00573 }
00574 }
00575
00576
00577
00578 for (long idx=noJentries-1; idx>=0; idx--)
00579 {
00580 bi = columnJentries[idx];
00581
00582 p_blockJ_pattern[bi] = 0;
00583
00584
00585 long Aij = columnJ.column_start_idx + idx*block_storage;
00586
00587
00588 Array::Copy(this->Columns_data,Aij,Atmp,0,block_storage);
00589
00590
00591 BlockArith->SubstSolveBlock(idd+bi*block_storage,cd+Aij);
00592
00593
00594
00595
00596 BlockArith->SubATBproduct(dd+Djj,atmp,cd+Aij);
00597 }
00598
00599
00600
00601 long maxJentries = columnJ.Entries-1;
00602 if (maxJentries<2)
00603 maxJentries = 2;
00604 if (columnJ.Entries>maxJentries)
00605 columnJ.Entries = maxJentries;
00606
00607
00608
00609
00610
00611
00612
00613
00614
00615
00616
00617
00618
00619
00620
00621
00622
00623
00624
00625
00626
00627
00628 }
00629
00630
00631 BlockArith->FactorizeBlock(dd+Djj);
00632
00633 if ((bj % newdotcnount) == 0)
00634 Write(".");
00635
00636 Djj += block_storage;
00637 eMT->act_block += block_size;
00638 if(eMT->break_flag)
00639 break;
00640 }
00641 }
00642 delete [] p_blockJ_pattern;
00643 delete [] Atmp;
00644 ComputeBlocks();
00645 }
00646
00647
00648
00649
00650 void SparseGridMtxLDL::SchurComplementFactorization(int fixed_blocks)
00651 {
00652 double* Atmp = new double[block_storage];
00653 long bi,bj,min_bi_J = 0;
00654
00655 long blocks_to_factor = n_blocks - fixed_blocks;
00656
00657
00658 long* p_blockJ_pattern = new long[n_blocks+1];
00659 memset(p_blockJ_pattern,0,(n_blocks+1)*sizeof(long));
00660
00661
00662
00663 BlockArith->zero_pivots = 0;
00664
00665 double* cd = this->Columns_data;
00666 double* atmp = Atmp;
00667 {
00668 double* dd = cd;
00669 double* idd = cd;
00670 long Djj = 0;
00671 for (bj=0; bj<n_blocks;bj++)
00672 {
00673 SparseGridColumn &columnJ = *Columns[bj];
00674 long noJentries = columnJ.Entries;
00675
00676
00677 long* columnJentries = columnJ.IndexesUfa->Items;
00678 double* pAkj = cd+columnJ.column_start_idx;
00679 double* pAij = pAkj;
00680
00681
00682 min_bi_J=bj;
00683
00684 if (bj<blocks_to_factor)
00685 {
00686 for (long i=noJentries-1; i>=0; i--)
00687 p_blockJ_pattern[min_bi_J=columnJentries[i]] = ~(i*block_storage);
00688 }
00689 else
00690 {
00691 for (long i=noJentries-1; i>=0; i--)
00692 if (columnJentries[i]<blocks_to_factor)
00693 p_blockJ_pattern[min_bi_J=columnJentries[i]] = ~(i*block_storage);
00694 }
00695
00696
00697 for (long idx_J=1; idx_J<noJentries; idx_J++)
00698 {
00699 pAij += block_storage;
00700 bi = columnJentries[idx_J];
00701
00702 SparseGridColumn &columnI = *Columns[bi];
00703 long noIentries = columnI.Entries;
00704
00705 if (noIentries>0)
00706 {
00707 double* pAki = cd+columnI.column_start_idx + (noIentries-1)*block_storage;
00708 long* columnIentries = columnI.IndexesUfa->Items;
00709 for (long* columnIentry = columnIentries+noIentries-1; columnIentry>=columnIentries; pAki-=block_storage)
00710 {
00711 long idx_K = p_blockJ_pattern[*columnIentry--];
00712 if (idx_K==0) continue;
00713 BlockArith->SubATBproduct(pAij,pAki,pAkj+ ~idx_K );
00714 }
00715 }
00716 }
00717
00718
00719
00720 for (long idx=0; idx<noJentries; idx++)
00721 {
00722 bi = columnJentries[idx];
00723
00724 p_blockJ_pattern[bi] = 0;
00725
00726 if (bi>=blocks_to_factor)
00727 break;
00728
00729
00730 long Aij = columnJ.column_start_idx + idx*block_storage;
00731
00732
00733 Array::Copy(this->Columns_data,Aij,Atmp,0,block_storage);
00734
00735
00736 BlockArith->SubstSolveBlock(idd+bi*block_storage,cd+Aij);
00737
00738
00739
00740
00741 BlockArith->SubATBproduct(dd+Djj,atmp,cd+Aij);
00742 }
00743
00744 if (bj<blocks_to_factor)
00745 BlockArith->FactorizeBlock(dd+Djj);
00746
00747
00748
00749
00750 Djj += block_storage;
00751 }
00752 }
00753
00754 delete [] p_blockJ_pattern;
00755 delete [] Atmp;
00756 ComputeBlocks();
00757 }
00758
00759
00760
00761 void SparseGridMtxLDL::Solve(double* b, double* x)
00762 {
00763 if (x!=b)
00764 Array::Copy(b,x,n);
00765 SolveLDL(x);
00766 }
00767
00768 void SparseGridMtxLDL::SolveLV(const LargeVector& f, LargeVector& r)
00769 {
00770 SolveLDL_node_perm(f, r);
00771 }
00772
00773 void SparseGridMtxLDL::SolveLDL_node_perm(const LargeVector& b, LargeVector& x)
00774 {
00775 if (node_order!=NULL)
00776 {
00777 b.GetPermuted_1Vector(tmp_vector_BS_nodes,node_order->perm->Items);
00778 SolveLDL_block_perm(*tmp_vector_BS_nodes,*tmp_vector_BS_nodes);
00779 tmp_vector_BS_nodes->GetPermutedVector(&x,node_order->perm->Items);
00780 }
00781 else
00782 {
00783 SolveLDL_block_perm(b,x);
00784 }
00785 }
00786
00787 void SparseGridMtxLDL::SolveLDL_block_perm(const LargeVector& b, LargeVector& x)
00788 {
00789 Solve(x.DataPtr(),b.DataPtr());
00790 }
00791
00792
00793 void SparseGridMtxLDL::SolveLDL(double* x,long fixed_blocks)
00794 {
00795 ForwardSubstL(x,fixed_blocks);
00796 SolveD(x,fixed_blocks);
00797 BackSubstLT(x,fixed_blocks);
00798 }
00799
00800
00801 void SparseGridMtxLDL::SubMultL12T(double* px, double* py,long fixed_blocks)
00802 {
00803 long schur_shift = n_blocks-fixed_blocks;
00804 if (fixed_blocks==0)
00805 return;
00806 long* ord = block_order->order->Items;
00807 double* cd = this->Columns_data;
00808 for (long bj=schur_shift; bj<n_blocks; bj++)
00809 {
00810 SparseGridColumn& columnJ = *Columns[bj];
00811
00812 long no = columnJ.Entries;
00813 if (no==0) continue;
00814
00815 double* dst = py+block_size*ord[bj];
00816 double* Aij = cd + columnJ.column_start_idx;
00817
00818 long* idxs = columnJ.IndexesUfa->Items;
00819
00820 for (long idx = 0; idx <no && *idxs<schur_shift; idx++,Aij += block_storage)
00821 BlockArith->SubMultTBlockByVector(Aij,px+block_size*ord[*(idxs++)],dst);
00822 }
00823 }
00824
00825
00826 void SparseGridMtxLDL::SubMultL12(double* px, double* py,long fixed_blocks)
00827 {
00828 long schur_shift = n_blocks-fixed_blocks;
00829 if (fixed_blocks==0)
00830 return;
00831 long* ord = block_order->order->Items;
00832 double* cd = this->Columns_data;
00833 for (long bj=schur_shift; bj<n_blocks; bj++)
00834 {
00835 SparseGridColumn& columnJ = *Columns[bj];
00836
00837 long no = columnJ.Entries;
00838 if (no==0) continue;
00839
00840 double* src = px+block_size*ord[bj];
00841 double* Aij = cd + columnJ.column_start_idx;
00842
00843 long* idxs = columnJ.IndexesUfa->Items;
00844
00845 for (long idx = 0; idx <no && *idxs<schur_shift; idx++,Aij += block_storage)
00846 BlockArith->SubMultBlockByVector(Aij,src,py+block_size*ord[*(idxs++)]);
00847 }
00848 }
00849
00850
00851 void SparseGridMtxLDL::ForwardSubstL(double* x,long fixed_blocks)
00852 {
00853 if (this->N()==0) return;
00854 long blocks_to_factor = n_blocks - fixed_blocks;
00855 long* ord = block_order->order->Items;
00856
00857 for (long bi=1; bi<blocks_to_factor; bi++)
00858 {
00859 SparseGridColumn& rowI = *Columns[bi];
00860 long no = rowI.Entries;
00861 if (no==0) continue;
00862 double* dst = x+block_size*ord[bi];
00863 double* Aij = Columns_data + rowI.column_start_idx;
00864
00865 long* idxs = rowI.IndexesUfa->Items;
00866
00867 for (long idx = 0; idx <no; idx++,Aij += block_storage)
00868 BlockArith->SubMultTBlockByVector(Aij,x+block_size*ord[*(idxs++)],dst);
00869 }
00870 }
00871
00872 void SparseGridMtxLDL::SolveD(double* x,long fixed_blocks)
00873 {
00874 if (this->N()==0) return;
00875 long* ord = block_order->order->Items;
00876
00877 double* Dii = Columns_data;
00878 long blocks_to_factor = n_blocks - fixed_blocks;
00879 for (long bi=0; bi<blocks_to_factor; bi++,Dii+=block_storage)
00880 BlockArith->SubstSolve(Dii,x+block_size*ord[bi]);
00881
00882 }
00883
00884 void SparseGridMtxLDL::BackSubstLT(double* x,long fixed_blocks)
00885 {
00886 if (this->N()==0) return;
00887 long* ord = block_order->order->Items;
00888 double* cd = this->Columns_data;
00889
00890 for (long bi=n_blocks-fixed_blocks-1; bi>=0; bi--)
00891 {
00892 SparseGridColumn& columnI = *Columns[bi];
00893 long no = columnI.Entries;
00894 if (no==0) continue;
00895 double* src = x+block_size*ord[bi];
00896 double* Aij = cd + columnI.column_start_idx;
00897 long* idxs = columnI.IndexesUfa->Items;
00898
00899 for (long idx = 0; idx <no; idx++,Aij += block_storage)
00900 BlockArith->SubMultBlockByVector(Aij,src,x+block_size*ord[*(idxs++)]);
00901 }
00902 }
00903
00904 void SparseGridMtxLDL::SolveA11(double* x,long fixed_blocks)
00905 {
00906 SolveLDL(x,fixed_blocks);
00907 }
00908
00909 void SparseGridMtxLDL::Sub_A21_A11inv(double* x,long fixed_blocks)
00910 {
00911 ForwardSubstL(x,fixed_blocks);
00912 SubMultL12T(x,x,fixed_blocks);
00913 SolveD(x,fixed_blocks);
00914 }
00915
00916 void SparseGridMtxLDL::Sub_A11inv_A12(double* x,long fixed_blocks)
00917 {
00918 SubMultL12(x,x,fixed_blocks);
00919 BackSubstLT(x,fixed_blocks);
00920 }
00921
00922 void SparseGridMtxLDL::WriteCondensedMatrixA22(double* a,Ordering* mcn,IntArrayList* lncn)
00923 {
00924 long blockSize = BlockSize();
00925 long i,j;
00926 long nbn = lncn->Count;
00927 long aux_bi_idx=0,aux_bj_idx=0;
00928 for (j=0; j<lncn->Count; j++)
00929 {
00930 long nj = mcn->perm->Items[lncn->Items[j]];
00931 long nbj = nj / blockSize;
00932 long nsj = nj % blockSize;
00933 long obj = block_order->perm->Items[nbj];
00934
00935 for (i=0; i<=j; i++)
00936 {
00937 long ni = mcn->perm->Items[lncn->Items[i]];
00938 long nbi = ni / blockSize;
00939 long nsi = ni % blockSize;
00940 long obi = block_order->perm->Items[nbi];
00941
00942 double val = GetValue(obi,obj,nsi,nsj,aux_bi_idx,aux_bj_idx);
00943 a[(j) + (i)*nbn] = val;
00944 a[(i) + (j)*nbn] = val;
00945 }
00946 }
00947 }
00948
00949 DSS_NAMESPASE_END