00001
00002
00003 #include "SparseConectivityMtx.h"
00004 #include "BiSection.h"
00005
00006 DSS_NAMESPASE_BEGIN
00007
00008 SparseConectivityMtxII::SparseConectivityMtxII(const SparseMatrixF& sm,char block_size)
00009 {
00010 n = (long)sm.neq/block_size + ( ((sm.neq % block_size) == 0) ? 0 : 1 );
00011 nonzeros = 0;
00012 ColumnsIndexes = new IntArrayList*[n];
00013 for (long bj=0; bj<n; bj++)
00014 ColumnsIndexes[bj] = new IntArrayList();
00015
00016 long* pmask = new long[n];
00017 memset(pmask,0,n*sizeof(long));
00018 IntArrayList* this_columnJ = NULL;
00019
00020
00021 for (ULONG j=0; j<sm.neq; j++)
00022 {
00023 long bj = j / block_size;
00024
00025 if (j % block_size == 0)
00026 this_columnJ = ColumnsIndexes[bj];
00027
00028 long lastbi=-1;
00029 for (ULONG ad = sm.Adr(j); ad<sm.Adr(j+1); ad++)
00030 {
00031 long i = sm.Ci(ad);
00032 long bi = i / block_size;
00033 if (bi==bj || bi==lastbi) continue;
00034 lastbi=bi;
00035
00036 if (pmask[bi]==0)
00037 {
00038 nonzeros++;
00039 this_columnJ->SortInsert(bi);
00040 pmask[bi] = 1;
00041 }
00042 }
00043
00044 if ((long)(j % block_size) == (long)(block_size-1))
00045 this_columnJ->SetIndexesTo(pmask,0);
00046 }
00047
00048
00049 for (ULONG i=0; i<sm.neq; i++)
00050 {
00051 long bi = i / block_size;
00052
00053 long lastbj=-1;
00054 for (ULONG ad = sm.Adr(i); ad<sm.Adr(i+1); ad++)
00055 {
00056 long j = sm.Ci(ad);
00057 long bj = j / block_size;
00058 if (bi==bj || bj==lastbj) continue;
00059 lastbj=bj;
00060
00061 if (ColumnsIndexes[bj]->AddIfNotLast(bi)>=0)
00062 nonzeros++;
00063 }
00064 }
00065
00066 if (pmask) delete [] pmask;
00067 pmask = NULL;
00068
00069 }
00070
00071 SparseConectivityMtxII::SparseConectivityMtxII(const SparseMatrixF& sm,Ordering* node_order,char block_size)
00072 {
00073 long* nprm = NULL;
00074 if (node_order)
00075 {
00076 nprm = node_order->perm->Items;
00077 n = node_order->order->Count/block_size;
00078 }
00079 else
00080 {
00081 n = (long)sm.neq/block_size + ( ((sm.neq % block_size) == 0) ? 0 : 1 );
00082 }
00083 nonzeros = 0;
00084 ColumnsIndexes = new IntArrayList*[n];
00085 for (long bj=0; bj<n; bj++)
00086 ColumnsIndexes[bj] = new IntArrayList();
00087
00088 long* pmask = new long[n];
00089 memset(pmask,0,n*sizeof(long));
00090 IntArrayList* this_columnJ = NULL;
00091
00092
00093 long old_bj = -1;
00094
00095
00096 for (ULONG j=0; j<sm.neq; j++)
00097 {
00098 long nj = node_order?node_order->perm->Items[j]:j;
00099 long bj = nj / block_size;
00100
00101 if (bj != old_bj)
00102 {
00103 if (this_columnJ) this_columnJ->SetIndexesTo(pmask,0);
00104 this_columnJ = ColumnsIndexes[bj];
00105 old_bj = bj;
00106 }
00107
00108 long lastbi=-1;
00109 for (ULONG ad = sm.Adr(j); ad<sm.Adr(j+1); ad++)
00110 {
00111 long i = sm.Ci(ad);
00112
00113 long ni = nprm?nprm[i]:i;
00114
00115 long bi = ni / block_size;
00116 if (bi==bj || bi==lastbi) continue;
00117 lastbi=bi;
00118
00119 if (pmask[bi]==0)
00120 {
00121 nonzeros++;
00122 this_columnJ->SortInsert(bi);
00123 pmask[bi] = 1;
00124 }
00125 }
00126 }
00127
00128
00129 for (ULONG i=0; i<sm.neq; i++)
00130 {
00131 long ni = nprm?nprm[i]:i;
00132
00133 long bi = ni / block_size;
00134
00135 long lastbj=-1;
00136 for (ULONG ad = sm.Adr(i); ad<sm.Adr(i+1); ad++)
00137 {
00138 long j = sm.Ci(ad);
00139 long nj = nprm?nprm[j]:j;
00140
00141 long bj = nj / block_size;
00142 if (bi==bj || bj==lastbj) continue;
00143 lastbj=bj;
00144
00145
00146
00147
00148 if (ColumnsIndexes[bj]->SortInsertIfNot(bi)>=0)
00149 nonzeros++;
00150 }
00151 }
00152
00153 if (pmask) delete [] pmask;
00154 pmask = NULL;
00155 }
00156
00157 SparseConectivityMtxII::SparseConectivityMtxII(SparseConectivityMtxII& mtx,Ordering* order)
00158 {
00159 n = mtx.N();
00160 ColumnsIndexes = new IntArrayList* [n];
00161 nonzeros = mtx.Nonzeros();
00162 long nj,idx,i;
00163 long* lengths = new long[n];
00164 memset(lengths,0,n*sizeof(long));
00165
00166
00167 for (nj=0; nj<n; nj++)
00168 {
00169 long j = order->order->Items[nj];
00170 IntArrayList& columnJ = *mtx.ColumnsIndexes[j];
00171
00172 for (idx = columnJ.Count-1; idx>=0; idx--)
00173 {
00174 long i = columnJ.Items[idx];
00175 long ni = order->perm->Items[i];
00176 if (ni>nj)
00177 lengths[ni]++;
00178 }
00179 }
00180
00181
00182 for (i=0; i<n; i++)
00183 {
00184 ColumnsIndexes[i] = new IntArrayList(lengths[i]);
00185 ColumnsIndexes[i]->Count = ColumnsIndexes[i]->Capacity;
00186 lengths[i] = 0;
00187 }
00188
00189
00190 for (nj=0; nj<n; nj++)
00191 {
00192 long j = order->order->Items[nj];
00193 IntArrayList& columnJ = *mtx.ColumnsIndexes[j];
00194
00195 long cnt = columnJ.Count;
00196 for (long idx = 0; idx<cnt; idx++)
00197 {
00198 long i = columnJ.Items[idx];
00199 long ni = order->perm->Items[i];
00200 if (ni>nj)
00201 ColumnsIndexes[ni]->Items[lengths[ni]++] = nj;
00202 }
00203 }
00204
00205 delete [] lengths;
00206 }
00207
00208
00209 SparseConectivityMtxII::SparseConectivityMtxII(const SparseConectivityMtxII& mtx)
00210 {
00211 n = mtx.N();
00212 ColumnsIndexes = new IntArrayList* [n];
00213 nonzeros = mtx.Nonzeros();
00214
00215
00216 for (long i=0; i<n; i++)
00217 ColumnsIndexes[i] = new IntArrayList(*mtx.ColumnsIndexes[i]);
00218 }
00219
00220 SparseConectivityMtxII::~SparseConectivityMtxII()
00221 {
00222 for (long i=0; i<n; i++)
00223 if (ColumnsIndexes[i]) {delete ColumnsIndexes[i]; ColumnsIndexes[i] = NULL;}
00224 delete [] ColumnsIndexes;
00225 ColumnsIndexes = NULL;
00226 }
00227
00228
00229 void SparseConectivityMtxII::AddGrow(long i, long j)
00230 {
00231 if (i<0 || j<0 || i>=n)
00232
00233 return;
00234
00235 IntArrayList* ColumnIndexes = ColumnsIndexes[i];
00236
00237 long myIndex=ColumnIndexes->BinarySearch(j);
00238 if ( myIndex < 0 )
00239 {
00240 ColumnIndexes->Insert(~myIndex,j);
00241 nonzeros++;
00242 }
00243 }
00244
00245
00246 long SparseConectivityMtxII::ColumnLength(long i)
00247 {
00248 return ColumnsIndexes[i]->Count;
00249 }
00250
00251 void SparseConectivityMtxII::GetCmpRows(long* &adr, long* &ci)
00252 {
00253 int n = N();
00254 adr = new long[n+1];
00255
00256 int i,length = 0;
00257 adr[0] = 0;
00258 for (i=0; i<n; i++)
00259 {
00260 IntArrayList* plist = GetIndexesAboveDiagonalInColumn(i);
00261 length += plist->Count;
00262 adr[i+1] = length;
00263 }
00264
00265 ci = new long[length];
00266 for (i=0; i<n; i++)
00267 {
00268 IntArrayList* plist = GetIndexesAboveDiagonalInColumn(i);
00269 Array::Copy(plist->Items,ci+adr[i],plist->Count);
00270 }
00271 }
00272
00273
00274 IntArrayList* SparseConectivityMtxII::GetOrder_Cuthill_McKee2()
00275 {
00276
00277 IntArrayList* order = new IntArrayList(n);
00278 order->Alloc();
00279 long* p_order = order->Items;
00280
00281
00282 long* p_node_level = new long[n];
00283 memset(p_node_level,0,n*sizeof(long));
00284
00285 long last_level_start = 0;
00286 long next_level_start = 0;
00287
00288 long last_level_count = 0;
00289 long next_level_count = 0;
00290
00291 long l=1;
00292 long k=0;
00293
00294 do
00295 {
00296 if (last_level_count == 0)
00297 {
00298 long r=-1,cl,min = n+1;
00299 for (long i=0; i<n; i++)
00300 if ((p_node_level[i]==0) && (cl=ColumnLength(i))<min) {min = cl;r=i;}
00301
00302 last_level_start = k;
00303
00304 if (r==-1)
00305 {
00306
00307 break;
00308 }
00309 p_node_level[r] = l;
00310 p_order[k++] = r;last_level_count=1;
00311 }
00312 next_level_start = k;
00313 next_level_count = 0;
00314
00315 for (long ilr = last_level_start; ilr<next_level_start; ilr++)
00316 {
00317 long lr = p_order[ilr];
00318
00319 IntArrayList& al= *ColumnsIndexes[lr];
00320 for (long idx = al.Count-1; idx>=0; idx--)
00321 {
00322 long neig = al.Items[idx];
00323 if (p_node_level[neig]!=0) continue;
00324 p_node_level[neig] = l+1;
00325 p_order[k++] = neig;next_level_count++;
00326 }
00327 }
00328 l++;
00329 last_level_count = next_level_count;
00330 last_level_start = next_level_start;
00331 }while(k<n);
00332
00333 delete [] p_node_level;
00334 return order;
00335 }
00336
00337 Ordering* SparseConectivityMtxII::GenerateMD(IntArrayList* fixed)
00338 {
00339 MD_Qqraph amd(this);
00340 return new Ordering(amd.GenerateMD(false,fixed));
00341 }
00342
00343 Ordering* SparseConectivityMtxII::GenerateAMD(IntArrayList* fixed)
00344 {
00345 MD_Qqraph amd(this);
00346 return new Ordering(amd.GenerateMD(true,fixed));
00347 }
00348
00349 Ordering* SparseConectivityMtxII::GenerateAMD_AA(IntArrayList* fixed)
00350 {
00351 MD_Qqraph amd(this);
00352 amd.aggressive_absorbtion = true;
00353 return new Ordering(amd.GenerateMD(true,fixed));
00354 }
00355
00356 Ordering* SparseConectivityMtxII::Get_Cuthill_McKee()
00357 {
00358 return new Ordering(this->GetOrder_Cuthill_McKee2());
00359 }
00360
00361 Ordering* SparseConectivityMtxII::Get_Reverse_Cuthill_McKee()
00362 {
00363 IntArrayList* order = this->GetOrder_Cuthill_McKee2();
00364 Array::Reverse(order->Items,order->Count);
00365 return new Ordering(order);
00366 }
00367
00368 Ordering* SparseConectivityMtxII::Get_Unity()
00369 {
00370 IntArrayList* order = new IntArrayList(n);
00371 order->InitIdentity();
00372 return new Ordering(order);
00373 }
00374
00375 Ordering* SparseConectivityMtxII::Get_RecursiveBiSection()
00376 {
00377 IntArrayList* order = new IntArrayList(n);
00378 order->InitIdentity();
00379
00380 CBiSection(this).RecurBiSectOrder(order);
00381
00382
00383
00384 return new Ordering(order);
00385 }
00386
00387 #ifdef _LINK_METIS_
00388 extern "C" void METIS_EdgeND(int *, int *, int *, int *, int *, int *, int *);
00389 extern "C" void METIS_NodeND(int *, int *, int *, int *, int *, int *, int *);
00390 #endif
00391
00392 Ordering* SparseConectivityMtxII::Get_MetisDiSection()
00393 {
00394 #ifdef _LINK_METIS_
00395 IntArrayList* order = new IntArrayList(n);
00396 IntArrayList* perm = new IntArrayList(n);
00397
00398 int n = N();
00399 long* adr;
00400 long* ci;
00401
00402 GetCmpRows(adr,ci);
00403
00404 int numflag = 0;
00405
00406
00407
00408
00409
00410
00411
00412
00413
00414
00415
00416
00417
00418
00419
00420 int optionsEdgeND[]=
00421 {
00422 0,
00423 3,
00424 1,
00425 2,
00426 0,
00427 0,
00428 0,
00429 1
00430 };
00431
00432 METIS_NodeND(&n,(int*)adr,(int*)ci,&numflag,optionsEdgeND,(int*)order->Items,(int*)perm->Items);
00433
00434 delete [] adr;
00435 delete [] ci;
00436
00437 return new Ordering(perm,order);
00438 #else
00439 Writeln(" The program was not linked with METIS library, use _LINK_METIS_ directive!");
00440 return NULL;
00441 #endif
00442 }
00443
00444 #ifdef _LINK_COLAMD_
00445 #include "colamd\colamd.h"
00446 #endif
00447
00448 Ordering* SparseConectivityMtxII::Get_ColAMD()
00449 {
00450 #ifdef _LINK_COLAMD_
00451 IntArrayList* order = new IntArrayList(n+1);
00452 order->Count = n;
00453
00454 int n = N();
00455 long* adr;
00456 long* ci;
00457
00458 GetCmpRows(adr,ci);
00459
00460 int stats [COLAMD_STATS];
00461
00462 symamd
00463 (
00464 n,
00465 (int*)ci,
00466 (int*)adr,
00467 (int*)order->Items,
00468 NULL,
00469 stats,
00470 calloc,
00471 free
00472 ) ;
00473
00474
00475
00476 delete [] adr;
00477 delete [] ci;
00478
00479 return new Ordering(order);
00480 #else
00481 Writeln(" The program was not linked with COLAMD sources, use _LINK_COLAMD_ directive!");
00482 return NULL;
00483 #endif
00484 }
00485
00486 void SparseConectivityMtxII::GenerateFillInPresorted(Ordering* ord)
00487 {
00488 MD_Qqraph amd(this);
00489 amd.keep_sorted_order = true;
00490 amd.GenerateMD(false,ord->order);
00491 }
00492
00493 Ordering* SparseConectivityMtxII::GetOrdering(Ordering::Type ord)
00494 {
00495 if (ord == Ordering::None)
00496 return Get_Unity();
00497 else
00498 if (ord == Ordering::ReverseCuthillMcKee)
00499 return Get_Reverse_Cuthill_McKee();
00500 else
00501 if (ord == Ordering::CuthillMcKee)
00502 return Get_Cuthill_McKee();
00503 else
00504 if (ord == Ordering::MinimumDegree)
00505 return GenerateMD();
00506 else
00507 if (ord == Ordering::ApproxMinimumDegree)
00508 return GenerateAMD();
00509 else
00510 if (ord == Ordering::MetisND)
00511 return Get_MetisDiSection();
00512 else
00513 if (ord == Ordering::ColAMD)
00514 return Get_ColAMD();
00515 else
00516 return NULL;
00517 }
00518
00519 IntArrayList* SparseConectivityMtxII::GetIndexesAboveDiagonalInColumn(long j)
00520 {
00521 return ColumnsIndexes[j];
00522 }
00523
00524 IntArrayList* SparseConectivityMtxII::DetachIndexesAboveDiagonalInColumn(long j)
00525 {
00526 IntArrayList* detach = ColumnsIndexes[j];
00527 ColumnsIndexes[j] = NULL;
00528 return detach;
00529 }
00530
00531 Ordering* SparseConectivityMtxII::GetPermutationAndPattern(Ordering::Type ord,IntArrayList* fixed)
00532 {
00533 Ordering* order = NULL;
00534 if (ord == Ordering::MinimumDegree)
00535 {
00536
00537 Writeln(" ordering : MinimumDegree");
00538 Write ("Symbolic QG factorization : ");
00539 clock_t start = MT.ClockStart();
00540 order = GenerateMD(fixed);
00541 Write(MT.MeasureClock(start));
00542
00543
00544
00545
00546 order->cm = new SparseConectivityMtxII(*this,order);
00547
00548 }
00549 else if (ord == Ordering::ApproxMinimumDegree)
00550 {
00551
00552 Writeln(" ordering : ApproxMinimumDegree");
00553 Write ("Symbolic QG factorization : ");
00554 order = GenerateAMD(fixed);
00555
00556
00557
00558 order->cm = new SparseConectivityMtxII(*this,order);
00559
00560 }
00561 else if (ord == Ordering::ApproxMinimumDegreeAA)
00562 {
00563
00564 Writeln(" ordering : ApproxMinimumDegree (aggressive absorbtion)");
00565 Write ("Symbolic QG factorization : ");
00566 order = GenerateAMD(fixed);
00567 GenerateFillInPresorted(order);
00568 order->cm = new SparseConectivityMtxII(*this,order);
00569 }
00570 else if (ord == Ordering::ApproxMinimumDegreeIncomplete)
00571 {
00572
00573
00574 Writeln(" ordering : ApproxMinimumDegreeIncomplete");
00575 Write ("Symbolic QG factorization : ");
00576 SparseConectivityMtxII* mtxA = new SparseConectivityMtxII(*this);
00577 order = GenerateAMD(fixed);
00578
00579
00580
00581 order->cm = new SparseConectivityMtxII(*mtxA,order);
00582 delete mtxA;
00583
00584 }
00585 else if (ord == Ordering::None)
00586 {
00587
00588 Writeln(" ordering : None");
00589 Write ("Symbolic QG factorization : ");
00590 order = Get_Unity();
00591 GenerateFillInPresorted(order);
00592
00593
00594
00595
00596 order->cm = new SparseConectivityMtxII(*this,order);
00597
00598 }
00599 else if (ord == Ordering::ReverseCuthillMcKee)
00600 {
00601 Writeln(" ordering : Reverse Cuthill-McKee");
00602 Write ("Symbolic QG factorization : ");
00603 clock_t start = MT.ClockStart();
00604 order = Get_Reverse_Cuthill_McKee();
00605 Write(MT.MeasureClock(start));Write("...");
00606 GenerateFillInPresorted(order);
00607 order->cm = new SparseConectivityMtxII(*this,order);
00608 }
00609 else if (ord == Ordering::CuthillMcKee)
00610 {
00611 Writeln(" ordering : Cuthill-McKee");
00612 Write ("Symbolic QG factorization : ");
00613 clock_t start = MT.ClockStart();
00614 order = Get_Cuthill_McKee();
00615 Write(MT.MeasureClock(start));Write("...");
00616 GenerateFillInPresorted(order);
00617 order->cm = new SparseConectivityMtxII(*this,order);
00618 }
00619 else if (ord == Ordering::NestedGraphBisection)
00620 {
00621 Writeln(" ordering : RecursiveBiSection");
00622 Write ("Symbolic QG factorization : ");
00623 clock_t start = MT.ClockStart();
00624 order = Get_RecursiveBiSection();
00625 Write(MT.MeasureClock(start));Write("...");
00626 GenerateFillInPresorted(order);
00627 order->cm = new SparseConectivityMtxII(*this,order);
00628 }
00629 else if (ord == Ordering::MetisND)
00630 {
00631 Writeln(" ordering : Metis (G.Karypis, V.Kumar)");
00632 Write ("Graph ordering optimization : ");
00633 clock_t start = MT.ClockStart();
00634 order = Get_MetisDiSection();
00635 Write(MT.MeasureClock(start));Write("...");
00636 GenerateFillInPresorted(order);
00637 order->cm = new SparseConectivityMtxII(*this,order);
00638 }
00639 else if (ord == Ordering::ColAMD)
00640 {
00641 Writeln(" ordering : ColAMD (S.I.Larimore, T.A.Davis)");
00642 Write ("Graph ordering optimization : ");
00643 clock_t start = MT.ClockStart();
00644 order = Get_ColAMD();
00645 Write(MT.MeasureClock(start));Write("...");
00646 GenerateFillInPresorted(order);
00647 order->cm = new SparseConectivityMtxII(*this,order);
00648 }
00649 else
00650 {
00651 Writeln("!!!! Incorrect ordering selected !!!!");
00652
00653 }
00654 return order;
00655 }
00656
00657
00658
00659
00660
00661
00662 MD_Qqraph::MD_Qqraph(SparseConectivityMtxII* mtx)
00663 {
00664 degrees = NULL;
00665 amd_w = NULL;
00666 pfixed = NULL;
00667
00668 MinDegA=MinDegB=MinDegC=vi_Min = 0;
00669 A = E = I = L = NULL;
00670 pLIdx = pPIdx = NULL;
00671 elements = NULL;
00672 no_elements = 0;
00673
00674 approximate_degree = false;
00675 keep_sorted_order = false;
00676 aggressive_absorbtion = false;
00677 this->mtx = mtx;
00678 this->n = mtx->N();
00679 n_1 = n-1;
00680 }
00681
00682 MD_Qqraph::~MD_Qqraph()
00683 {
00684 if (degrees) {delete [] degrees;degrees=NULL;}
00685 if (elements) {delete elements;elements = NULL;}
00686 if (amd_w) {delete [] amd_w;amd_w = NULL;}
00687
00688 for (long i=0; i<n; i++)
00689 {
00690
00691 if (E[i]) {delete E[i];E[i] = NULL;}
00692 if (L[i]) {
00693 delete L[i];L[i] = NULL;
00694 }
00695 if (I[i]) {delete I[i];I[i] = NULL;}
00696 }
00697
00698 delete [] E; E = NULL;
00699 delete [] L; L = NULL;
00700 delete [] I; I = NULL;
00701 }
00702
00703 void MD_Qqraph::Insert(IntArrayList& al, long i)
00704 {
00705 long myIndex = al.BinarySearch(i);
00706 if (myIndex<0)
00707 al.Insert(~myIndex,i);
00708 }
00709
00710 void MD_Qqraph::ComputeAMD_w_Le_Lp(long i)
00711 {
00712 IntArrayList& el = *E[i];
00713 for (long ie=el.Count-1; ie>=0; ie--)
00714 {
00715 long e = el[ie];
00716 if (amd_w[e]<0) amd_w[e] = L[e]->Count;
00717 amd_w[e]-=I[i]->Count;
00718
00719 if (aggressive_absorbtion && amd_w[e]==0)
00720 {
00721 el.RemoveAt(ie);
00722 amd_w[e] = -1;
00723 }
00724 }
00725 }
00726
00727 void MD_Qqraph::ClearAMD_w(long i)
00728 {
00729 IntArrayList& el = *E[i];
00730 for (long ie=el.Count-1; ie>=0; ie--)
00731 amd_w[el[ie]] = -1;
00732 }
00733
00734 long MD_Qqraph::ApproximateMinumumDegree(long i,IntArrayList& Lp)
00735 {
00736 IntArrayList &Ali = *A[i];
00737 IntArrayList &Eli = *E[i];
00738
00739 long idx;
00740
00741 long ne = Eli.Count;
00742 long* pe = Eli.Items;
00743 {
00744 long Aic = 0;
00745
00746 Aic += Ali.CountWithoutMask(pPIdx);
00747
00748 long Lpc = 0;
00749
00750 for(idx=0; idx<Lp.Count; idx++){long j = Lp[idx]; if (pPIdx[j]==0) Lpc++;}
00751
00752 long d=min(n - this->no_elements,degrees[i]+Lpc);
00753
00754 long Les = 0;
00755
00756 for(idx=ne-1; idx>=0; idx--)
00757 {
00758 long e = pe[idx];
00759 if (amd_w[e]<0) Les += L[e]->Count;
00760 else Les += amd_w[e];
00761 }
00762
00763 d = min(d,Aic + Lpc + Les);
00764
00765 if (MinDegB>d && IsFree(i))
00766 {
00767 MinDegB = d;
00768 vi_Min = i;
00769 }
00770 return d;
00771 }
00772 }
00773
00774
00775
00776
00777 long MD_Qqraph::ExternalDegree(long i)
00778 {
00779 long idx,idx2;
00780 long d = 0;
00781 long na = A[i]->Count;
00782 long ne = E[i]->Count;
00783 long* pa = A[i]->Items,*pe = E[i]->Items;
00784 {
00785
00786 for (idx=na-1; idx>=0; idx--)
00787 {
00788 long a = pa[idx];
00789 if(pPIdx[a]==0)
00790 {
00791 pLIdx[a] = 1;
00792 d++;
00793 }
00794 }
00795
00796
00797 for (idx=ne-1; idx>=0; idx--)
00798 {
00799 long e = pe[idx];
00800
00801 IntArrayList& Lp = *L[e];
00802 for(idx2=0; idx2<Lp.Count; idx2++)
00803 {
00804 long j = Lp[idx2];
00805 if (pPIdx[j]==0 && pLIdx[j]!=1)
00806 {
00807 pLIdx[j] = 1;
00808 d++;
00809 }
00810 }
00811 }
00812
00813
00814 A[i]->SetIndexesTo(pLIdx,0);
00815
00816
00817 for (idx=ne-1; idx>=0; idx--)
00818 {
00819 long e = pe[idx];
00820
00821 IntArrayList& Lp = *L[e];
00822 for(idx2=0; idx2<Lp.Count; idx2++)
00823 {
00824 long j = Lp[idx2];
00825 if(pPIdx[j]==0) pLIdx[j] = 0;
00826 }
00827 }
00828
00829 if (MinDegB>d && IsFree(i))
00830 {
00831 MinDegB = d;
00832 vi_Min = i;
00833 }
00834 }
00835 return d;
00836 }
00837
00838
00839
00840 void MD_Qqraph::Eliminate(long p)
00841 {
00842 if (I[p] == NULL)
00843 return;
00844
00845 long ie,it,ii;
00846
00847
00848 IntArrayList &Ip = *I[p];
00849
00850 Ip.SetIndexesTo(pLIdx,1);
00851
00852
00853 IntArrayList& Lp = *A[p];A[p] = NULL;
00854 Lp.RemoveMarkByBattern(pLIdx);
00855
00856 IntArrayList &Ep = *E[p];
00857
00858 for(ie = Ep.Count-1; ie>=0; ie--)
00859 {
00860 long e = Ep.Items[ie];
00861 pLIdx[e] = 1;
00862
00863 IntArrayList& Le = *L[e];
00864 for(it=0; it<Le.Count; it++)
00865 {
00866 long i = Le[it];
00867 if (pLIdx[i]==0)
00868 {
00869 Lp.SortInsert(i);
00870 pLIdx[i]=1;
00871 }
00872 }
00873 }
00874
00875
00876 L[p] = &Lp;
00877
00878 for(it=0; it<Lp.Count; it++)
00879 {
00880 long i = Lp[it];
00881 if (I[i]==NULL) continue;
00882
00883 A[i]->RemoveByBattern(pLIdx);
00884
00885 E[i]->RemoveByBattern(pLIdx);
00886 Insert(*E[i],p);
00887
00888
00889 if (!keep_sorted_order && approximate_degree)
00890 ComputeAMD_w_Le_Lp(i);
00891 }
00892
00893
00894 Lp. SetIndexesTo(pLIdx,0);
00895 I[p]->SetIndexesTo(pLIdx,0);
00896 E[p]->SetIndexesTo(pLIdx,0);
00897
00898 if (!keep_sorted_order)
00899 {
00900
00901
00902 for(it=0; it<Lp.Count; it++)
00903 {
00904 long i = Lp[it];
00905 if (I[i] == NULL) continue;
00906 I[i]->SetIndexesTo(pPIdx,1);
00907 if (approximate_degree)
00908 this->degrees[i] = ApproximateMinumumDegree(i,Lp);
00909 else
00910 this->degrees[i] = ExternalDegree(i);
00911 I[i]->SetIndexesTo(pPIdx,0);
00912 }
00913
00914 if (approximate_degree)
00915
00916 for(it=0; it<Lp.Count; it++)
00917 {
00918 long i = Lp[it];
00919 if (I[i]!=NULL) ClearAMD_w(i);
00920 }
00921
00922
00923 SupervariablesDetection(Lp);
00924 }
00925
00926 IntArrayList* Lpfull = new IntArrayList(Lp);
00927
00928
00929
00930 for(ii = Ip.Count-1; ii>=0; ii--)
00931 {
00932 long k = Ip.Items[ii];
00933 Insert(*Lpfull,k);
00934 }
00935 Lpfull->TrimToSize();
00936 mtx->ColumnsIndexes[p] = Lpfull;
00937
00938 for(ii = Ip.Count-1; ii>=0; ii--)
00939 {
00940 long k = Ip.Items[ii];
00941 elements->Items[no_elements++] = k;
00942 if (k != p)
00943 mtx->ColumnsIndexes[k]->Fill(*Lpfull);
00944 }
00945
00946 variables.Remove(p);
00947
00948
00949 if (E[p]) {delete E[p];E[p] = NULL;}
00950 if (I[p]) {delete I[p];I[p] = NULL;}
00951 }
00952
00953 bool MD_Qqraph::IsIndistinguishable(long i,long j)
00954 {
00955 if (A[i]->Count!=A[j]->Count) return false;
00956 if (E[i]->Count!=E[j]->Count) return false;
00957
00958 bool ind = true;
00959
00960 pPIdx[j] = pPIdx[i] = 1;
00961 A[i]->SetIndexesTo(pPIdx,1);
00962 E[i]->SetIndexesTo(pPIdx,2);
00963
00964 A[j]->TestSetIndexesTo(pPIdx,1,0,ind);
00965 E[j]->TestSetIndexesTo(pPIdx,2,0,ind);
00966
00967 pPIdx[j] = pPIdx[i] = 0;
00968 if (ind)
00969 return true;
00970
00971 A[i]->SetIndexesTo(pPIdx,0);
00972 E[i]->SetIndexesTo(pPIdx,0);
00973 return false;
00974 }
00975
00976 void MD_Qqraph::SupervariablesDetection(IntArrayList& Lp)
00977 {
00978 if (keep_sorted_order)
00979 return;
00980
00981 long idx,ix;
00982
00983
00984
00985 for(idx=0; idx<Lp.Count; idx++)
00986 {
00987 long i = Lp[idx];
00988 if (I[i]==NULL || !IsFree(i)) continue;
00989 long key = Hash(i);
00990 ht.AddValue(key,i);
00991 }
00992
00993 for (idx = ht.occupied_buckets.Count-1; idx>=0; idx--)
00994 {
00995 IntArrayList& al = *ht.buckets[ht.occupied_buckets.Items[idx]];
00996 if (al.Count>=2)
00997 {
00998 hash_parents.Alloc(al.Count);
00999 for (ix=0; ix<al.Count; ix++)
01000 {
01001 long i=al[ix];
01002 if (hash_parents[ix]<0) continue;
01003 for (long jx=ix+1; jx<al.Count; jx++)
01004 {
01005 if (hash_parents[jx]<0) continue;
01006 long j=al[jx];
01007 if (IsIndistinguishable(i,j))
01008 hash_parents[jx] = ~i;
01009 }
01010 }
01011
01012 for (ix=0; ix<al.Count; ix++)
01013 {
01014 if (hash_parents[ix]<0)
01015 {
01016 long i=~hash_parents[ix];
01017 long j=al[ix];
01018
01019
01020 IntArrayList &Ij = *I[j];
01021 for(long ij = Ij.Count-1; ij>=0; ij--)
01022 I[i]->SortInsert(Ij.Items[ij]);
01023
01024 degrees[i] -= I[j]->Count;
01025 if (MinDegB>degrees[i])
01026 {
01027 MinDegB = degrees[i];
01028 vi_Min = i;
01029 }
01030
01031 variables.Remove(j);
01032 if (I[j]) {delete I[j];I[j] = NULL;}
01033
01034
01035 if (E[j]) {delete E[j];E[j] = NULL;}
01036 }
01037 }
01038 }
01039 }
01040 ht.Clear();
01041 }
01042
01043
01044
01045
01046
01047
01048
01049 IntArrayList* MD_Qqraph::GenerateMD(bool approximate_degree,IntArrayList* fixed )
01050 {
01051 this->approximate_degree = approximate_degree;
01052 no_elements = 0;
01053
01054 ht.Init(n);
01055
01056
01057
01058
01059 long* permutation = new long[n];
01060
01061
01062
01063 long* vlasts = new long[n];
01064 long* vnexts = new long[n];
01065
01066
01067 long no_fixed = 0;
01068 if (fixed)
01069 {
01070 pfixed = new long[n];
01071 Array::Clear(pfixed,0,n);
01072 fixed->SetIndexesTo(pfixed,1);
01073 no_fixed = fixed->Count;
01074 }
01075
01076 elements = new IntArrayList(n);
01077 elements->Alloc();
01078 long* pivot_pattern = new long[n];
01079 memset(pivot_pattern,0,n*sizeof(long));
01080 degrees = new long[n];
01081 amd_w = new long[n];
01082
01083 MinDegA = n;
01084 MinDegB = n;
01085 MinDegC = n;
01086
01087 this->E = new IntArrayList*[n];
01088 this->L = new IntArrayList*[n];
01089 this->A = mtx->ColumnsIndexes;
01090 this->I = new IntArrayList*[n];
01091
01092 for (long k=0; k<n; k++)
01093 {
01094 I[k] = new IntArrayList(4);I[k]->Add(k);
01095 E[k] = new IntArrayList(4);
01096
01097 L[k] = NULL;
01098
01099 vlasts[k] = k-1;
01100 vnexts[k] = k+1;
01101
01102 degrees[k] = A[k]->Count;
01103 if (degrees[k]>n)
01104 {
01105 mtx->Writeln("Connectivity matrix is corrupt !");
01106 return NULL;
01107 }
01108
01109 if (degrees[k]<MinDegA && IsFree(k))
01110 MinDegA = degrees[k];
01111 permutation[k] = 0;
01112 amd_w[k] = -1;
01113 }
01114 vnexts[n-1] = -1;
01115
01116
01117
01118
01119
01120
01121
01122 long* pLIdx = permutation,*pPIdx = pivot_pattern,*vl = vlasts, *vn = vnexts;
01123 {
01124 this->pLIdx = pLIdx;
01125 this->pPIdx = pPIdx;
01126 variables.Init(vl,vn,n);
01127
01128 if (keep_sorted_order)
01129 {
01130 for (long i=0;i<n;i++)
01131 Eliminate(fixed->Items[i]);
01132 }
01133 else
01134 {
01135 while (no_elements<n)
01136 {
01137 long MinimDeg = min(MinDegB,MinDegA);
01138 MinDegB = n;
01139 MinDegC = n;
01140
01141
01142 for (long vi=variables.first;vi>=0;)
01143 {
01144
01145 if (IsFree(vi) && (MinimDeg >= degrees[vi]))
01146 {
01147
01148 Eliminate(vi);
01149 if (MinDegB<MinimDeg)
01150 {
01151 MinimDeg = MinDegB;
01152 vi = vi_Min;
01153 continue;
01154 }
01155 }
01156 else
01157 {
01158 if (IsFree(vi) && (degrees[vi]<MinDegC))
01159 MinDegC = degrees[vi];
01160 }
01161 vi=variables.next[vi];
01162 }
01163 MinDegA = MinDegC;
01164
01165 if (no_fixed && (no_elements==n-no_fixed))
01166 {
01167 delete [] pfixed;pfixed = NULL;
01168 no_fixed = 0;
01169 }
01170
01171
01172
01173
01174
01175
01176 }
01177 }
01178 }
01179
01180 delete [] pfixed;pfixed = NULL;
01181
01182 delete [] permutation;
01183 delete [] vlasts;
01184 delete [] vnexts;
01185
01186 delete [] pivot_pattern;
01187
01188 IntArrayList* order = elements;
01189 elements = NULL;
01190
01191 return order;
01192 }
01193
01194
01195 DSS_NAMESPASE_END