00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061
00062
00063 #define spINSIDE_SPARSE
00064 #include "spConfig.h"
00065 #include "spmatrix.h"
00066 #include "spDefs.h"
00067 #include "spmalloc.h"
00068
00069 static InitializeElementBlocks();
00070 static RecordAllocation();
00071 static AllocateBlockOfAllocationList();
00072
00073
00074
00075
00076
00077
00078
00079
00080
00081
00082
00083
00084
00085
00086
00087
00088
00089
00090
00091
00092
00093
00094
00095
00096
00097
00098
00099
00100
00101
00102
00103
00104
00105
00106
00107
00108
00109
00110
00111
00112 char *
00113 spCreate( Size, Complex, pError )
00114
00115 int Size, *pError;
00116 SPBOOLEAN Complex;
00117 {
00118 register unsigned SizePlusOne;
00119 register MatrixPtr Matrix;
00120 register int I;
00121 int AllocatedSize;
00122
00123
00124
00125 *pError = spOKAY;
00126
00127
00128 if ((Size < 0) OR (Size == 0 AND NOT EXPANDABLE))
00129 { *pError = spPANIC;
00130 return NULL;
00131 }
00132
00133
00134 #if NOT spCOMPLEX
00135 if (Complex)
00136 { *pError = spPANIC;
00137 return NULL;
00138 }
00139 #endif
00140 #if NOT REAL
00141 if (NOT Complex)
00142 { *pError = spPANIC;
00143 return NULL;
00144 }
00145 #endif
00146
00147
00148 AllocatedSize = MAX( Size, MINIMUM_ALLOCATED_SIZE );
00149 SizePlusOne = (unsigned)(AllocatedSize + 1);
00150
00151 if ((Matrix = SPALLOC(struct MatrixFrame, 1)) == NULL)
00152 { *pError = spNO_MEMORY;
00153 return NULL;
00154 }
00155
00156
00157 Matrix->ID = SPARSE_ID;
00158 Matrix->Complex = Complex;
00159 Matrix->PreviousMatrixWasComplex = Complex;
00160 Matrix->Factored = NO;
00161 Matrix->Elements = 0;
00162 Matrix->Error = *pError;
00163 Matrix->Fillins = 0;
00164 Matrix->Reordered = NO;
00165 Matrix->NeedsOrdering = YES;
00166 Matrix->NumberOfInterchangesIsOdd = NO;
00167 Matrix->Partitioned = NO;
00168 Matrix->RowsLinked = NO;
00169 Matrix->InternalVectorsAllocated = NO;
00170 Matrix->SingularCol = 0;
00171 Matrix->SingularRow = 0;
00172 Matrix->Size = Size;
00173 Matrix->AllocatedSize = AllocatedSize;
00174 Matrix->ExtSize = Size;
00175 Matrix->AllocatedExtSize = AllocatedSize;
00176 Matrix->CurrentSize = 0;
00177 Matrix->ExtToIntColMap = NULL;
00178 Matrix->ExtToIntRowMap = NULL;
00179 Matrix->IntToExtColMap = NULL;
00180 Matrix->IntToExtRowMap = NULL;
00181 Matrix->MarkowitzRow = NULL;
00182 Matrix->MarkowitzCol = NULL;
00183 Matrix->MarkowitzProd = NULL;
00184 Matrix->DoCmplxDirect = NULL;
00185 Matrix->DoRealDirect = NULL;
00186 Matrix->Intermediate = NULL;
00187 Matrix->RelThreshold = DEFAULT_THRESHOLD;
00188 Matrix->AbsThreshold = 0.0;
00189
00190 Matrix->TopOfAllocationList = NULL;
00191 Matrix->RecordsRemaining = 0;
00192 Matrix->ElementsRemaining = 0;
00193 Matrix->FillinsRemaining = 0;
00194
00195 RecordAllocation( Matrix, (char *)Matrix );
00196 if (Matrix->Error == spNO_MEMORY) goto MemoryError;
00197
00198
00199 Matrix->TrashCan.Real = 0.0;
00200 #if spCOMPLEX
00201 Matrix->TrashCan.Imag = 0.0;
00202 #endif
00203 Matrix->TrashCan.Row = 0;
00204 Matrix->TrashCan.Col = 0;
00205 Matrix->TrashCan.NextInRow = NULL;
00206 Matrix->TrashCan.NextInCol = NULL;
00207 #if INITIALIZE
00208 Matrix->TrashCan.pInitInfo = NULL;
00209 #endif
00210
00211
00212 SPCALLOC( Matrix->Diag, ElementPtr, SizePlusOne);
00213 if (Matrix->Diag == NULL)
00214 goto MemoryError;
00215
00216
00217 SPCALLOC( Matrix->FirstInCol, ElementPtr, SizePlusOne);
00218 if (Matrix->FirstInCol == NULL)
00219 goto MemoryError;
00220
00221
00222 SPCALLOC( Matrix->FirstInRow, ElementPtr, SizePlusOne);
00223 if (Matrix->FirstInRow == NULL)
00224 goto MemoryError;
00225
00226
00227 if (( Matrix->IntToExtColMap = SPALLOC(int, SizePlusOne)) == NULL)
00228 goto MemoryError;
00229
00230
00231 if (( Matrix->IntToExtRowMap = SPALLOC(int, SizePlusOne)) == NULL)
00232 goto MemoryError;
00233
00234
00235 for (I = 1; I <= AllocatedSize; I++)
00236 { Matrix->IntToExtRowMap[I] = I;
00237 Matrix->IntToExtColMap[I] = I;
00238 }
00239
00240 #if TRANSLATE
00241
00242 if (( Matrix->ExtToIntColMap = SPALLOC(int, SizePlusOne)) == NULL)
00243 goto MemoryError;
00244
00245
00246 if (( Matrix->ExtToIntRowMap = SPALLOC(int, SizePlusOne)) == NULL)
00247 goto MemoryError;
00248
00249
00250 for (I = 1; I <= AllocatedSize; I++)
00251 { Matrix->ExtToIntColMap[I] = -1;
00252 Matrix->ExtToIntRowMap[I] = -1;
00253 }
00254 Matrix->ExtToIntColMap[0] = 0;
00255 Matrix->ExtToIntRowMap[0] = 0;
00256 #endif
00257
00258
00259 InitializeElementBlocks( Matrix, SPACE_FOR_ELEMENTS*AllocatedSize,
00260 SPACE_FOR_FILL_INS*AllocatedSize );
00261 if (Matrix->Error == spNO_MEMORY)
00262 goto MemoryError;
00263
00264 return (char *)Matrix;
00265
00266 MemoryError:
00267
00268
00269
00270 *pError = spNO_MEMORY;
00271 spDestroy( (char *)Matrix);
00272 return NULL;
00273 }
00274
00275
00276
00277
00278
00279
00280
00281
00282
00283
00284
00285
00286
00287
00288
00289
00290
00291
00292
00293
00294
00295
00296
00297
00298
00299
00300
00301
00302
00303
00304
00305
00306 ElementPtr
00307 spcGetElement( Matrix )
00308
00309 MatrixPtr Matrix;
00310 {
00311 ElementPtr pElement;
00312
00313
00314
00315
00316 if (Matrix->ElementsRemaining == 0)
00317 { pElement = SPALLOC(struct MatrixElement, ELEMENTS_PER_ALLOCATION);
00318 RecordAllocation( Matrix, (char *)pElement );
00319 if (Matrix->Error == spNO_MEMORY) return NULL;
00320 Matrix->ElementsRemaining = ELEMENTS_PER_ALLOCATION;
00321 Matrix->NextAvailElement = pElement;
00322 }
00323
00324
00325 Matrix->ElementsRemaining--;
00326 return Matrix->NextAvailElement++;
00327
00328 }
00329
00330
00331
00332
00333
00334
00335
00336
00337
00338
00339
00340
00341
00342
00343
00344
00345
00346
00347
00348
00349
00350
00351
00352
00353
00354
00355
00356
00357
00358
00359
00360
00361
00362
00363
00364
00365
00366
00367 static
00368 InitializeElementBlocks( Matrix, InitialNumberOfElements,
00369 NumberOfFillinsExpected )
00370
00371 MatrixPtr Matrix;
00372 int InitialNumberOfElements, NumberOfFillinsExpected;
00373 {
00374 ElementPtr pElement;
00375
00376
00377
00378
00379 pElement = SPALLOC(struct MatrixElement, InitialNumberOfElements);
00380 RecordAllocation( Matrix, (char *)pElement );
00381 if (Matrix->Error == spNO_MEMORY) return 0;
00382 Matrix->ElementsRemaining = InitialNumberOfElements;
00383 Matrix->NextAvailElement = pElement;
00384
00385
00386 pElement = SPALLOC(struct MatrixElement, NumberOfFillinsExpected);
00387 RecordAllocation( Matrix, (char *)pElement );
00388 if (Matrix->Error == spNO_MEMORY) return 0;
00389 Matrix->FillinsRemaining = NumberOfFillinsExpected;
00390 Matrix->NextAvailFillin = pElement;
00391
00392
00393 Matrix->FirstFillinListNode = SPALLOC(struct FillinListNodeStruct,1);
00394 RecordAllocation( Matrix, (char *)Matrix->FirstFillinListNode );
00395 if (Matrix->Error == spNO_MEMORY) return 0;
00396 Matrix->LastFillinListNode = Matrix->FirstFillinListNode;
00397
00398 Matrix->FirstFillinListNode->pFillinList = pElement;
00399 Matrix->FirstFillinListNode->NumberOfFillinsInList =NumberOfFillinsExpected;
00400 Matrix->FirstFillinListNode->Next = NULL;
00401
00402 return 0;
00403 }
00404
00405
00406
00407
00408
00409
00410
00411
00412
00413
00414
00415
00416
00417
00418
00419
00420
00421
00422
00423
00424
00425
00426
00427
00428
00429
00430
00431
00432
00433 ElementPtr
00434 spcGetFillin( Matrix )
00435
00436 MatrixPtr Matrix;
00437 {
00438 struct FillinListNodeStruct *pListNode;
00439 ElementPtr pFillins;
00440
00441
00442
00443 #if NOT STRIP OR LINT
00444 if (Matrix->FillinsRemaining == 0)
00445 return spcGetElement( Matrix );
00446 #endif
00447 #if STRIP OR LINT
00448
00449 if (Matrix->FillinsRemaining == 0)
00450 { pListNode = Matrix->LastFillinListNode;
00451
00452
00453 if (pListNode->Next != NULL)
00454 { Matrix->LastFillinListNode = pListNode = pListNode->Next;
00455 Matrix->FillinsRemaining = pListNode->NumberOfFillinsInList;
00456 Matrix->NextAvailFillin = pListNode->pFillinList;
00457 }
00458 else
00459 {
00460
00461 pFillins = SPALLOC(struct MatrixElement, ELEMENTS_PER_ALLOCATION);
00462 RecordAllocation( Matrix, (char *)pFillins );
00463 if (Matrix->Error == spNO_MEMORY) return NULL;
00464 Matrix->FillinsRemaining = ELEMENTS_PER_ALLOCATION;
00465 Matrix->NextAvailFillin = pFillins;
00466
00467
00468 pListNode->Next = SPALLOC(struct FillinListNodeStruct,1);
00469 RecordAllocation( Matrix, (char *)pListNode->Next );
00470 if (Matrix->Error == spNO_MEMORY) return NULL;
00471 Matrix->LastFillinListNode = pListNode = pListNode->Next;
00472
00473 pListNode->pFillinList = pFillins;
00474 pListNode->NumberOfFillinsInList = ELEMENTS_PER_ALLOCATION;
00475 pListNode->Next = NULL;
00476 }
00477 }
00478 #endif
00479
00480
00481 Matrix->FillinsRemaining--;
00482 return Matrix->NextAvailFillin++;
00483 }
00484
00485
00486
00487
00488
00489
00490
00491
00492
00493
00494
00495
00496
00497
00498
00499
00500
00501
00502
00503
00504
00505
00506
00507
00508
00509
00510 static
00511 RecordAllocation( Matrix, AllocatedPtr )
00512
00513 MatrixPtr Matrix;
00514 char *AllocatedPtr;
00515 {
00516
00517
00518
00519
00520
00521 if (AllocatedPtr == NULL)
00522 { Matrix->Error = spNO_MEMORY;
00523 return 0;
00524 }
00525
00526
00527 if (Matrix->RecordsRemaining == 0)
00528 { AllocateBlockOfAllocationList( Matrix );
00529 if (Matrix->Error == spNO_MEMORY)
00530 { SPFREE(AllocatedPtr);
00531 return 0;
00532 }
00533 }
00534
00535
00536 (++Matrix->TopOfAllocationList)->AllocatedPtr = AllocatedPtr;
00537 Matrix->RecordsRemaining--;
00538 return 0;
00539
00540 }
00541
00542
00543
00544
00545
00546
00547
00548
00549
00550
00551
00552
00553
00554
00555
00556
00557
00558
00559
00560
00561
00562
00563
00564
00565
00566
00567 static
00568 AllocateBlockOfAllocationList( Matrix )
00569
00570 MatrixPtr Matrix;
00571 {
00572 register int I;
00573 register AllocationListPtr ListPtr;
00574
00575
00576
00577 ListPtr = SPALLOC(struct AllocationRecord, (ELEMENTS_PER_ALLOCATION+1));
00578 if (ListPtr == NULL)
00579 { Matrix->Error = spNO_MEMORY;
00580 return 0;
00581 }
00582
00583
00584
00585
00586 ListPtr->NextRecord = Matrix->TopOfAllocationList;
00587 Matrix->TopOfAllocationList = ListPtr;
00588 ListPtr += ELEMENTS_PER_ALLOCATION;
00589 for (I = ELEMENTS_PER_ALLOCATION; I > 0; I--)
00590 { ListPtr->NextRecord = ListPtr - 1;
00591 ListPtr--;
00592 }
00593
00594
00595 Matrix->TopOfAllocationList->AllocatedPtr = (char *)ListPtr;
00596 Matrix->RecordsRemaining = ELEMENTS_PER_ALLOCATION;
00597
00598 return 0;
00599 }
00600
00601
00602
00603
00604
00605
00606
00607
00608
00609
00610
00611
00612
00613
00614
00615
00616
00617
00618
00619
00620
00621
00622
00623
00624
00625
00626
00627
00628 void
00629 spDestroy( eMatrix )
00630
00631 register char *eMatrix;
00632 {
00633 MatrixPtr Matrix = (MatrixPtr)eMatrix;
00634 register AllocationListPtr ListPtr, NextListPtr;
00635
00636
00637
00638 ASSERT( IS_SPARSE( Matrix ) );
00639
00640
00641 SPFREE( Matrix->IntToExtColMap );
00642 SPFREE( Matrix->IntToExtRowMap );
00643 SPFREE( Matrix->ExtToIntColMap );
00644 SPFREE( Matrix->ExtToIntRowMap );
00645 SPFREE( Matrix->Diag );
00646 SPFREE( Matrix->FirstInRow );
00647 SPFREE( Matrix->FirstInCol );
00648 SPFREE( Matrix->MarkowitzRow );
00649 SPFREE( Matrix->MarkowitzCol );
00650 SPFREE( Matrix->MarkowitzProd );
00651 SPFREE( Matrix->DoCmplxDirect );
00652 SPFREE( Matrix->DoRealDirect );
00653 SPFREE( Matrix->Intermediate );
00654
00655
00656
00657
00658
00659 ListPtr = Matrix->TopOfAllocationList;
00660 while (ListPtr != NULL )
00661 {
00662 char *LocPtr;
00663
00664
00665
00666
00667
00668
00669
00670
00671
00672 NextListPtr = ListPtr->NextRecord;
00673
00674 LocPtr=ListPtr->AllocatedPtr;
00675 ListPtr->AllocatedPtr= NULL;
00676 if ( LocPtr != NULL) free(LocPtr);
00677 ListPtr = NextListPtr;
00678 }
00679 return;
00680 }
00681
00682
00683
00684
00685
00686
00687
00688
00689
00690
00691
00692
00693
00694
00695
00696
00697
00698
00699
00700
00701 int
00702 spError( eMatrix )
00703
00704 char *eMatrix;
00705 {
00706
00707
00708 if (eMatrix != NULL)
00709 { ASSERT(((MatrixPtr)eMatrix)->ID == SPARSE_ID);
00710 return ((MatrixPtr)eMatrix)->Error;
00711 }
00712 else return spNO_MEMORY;
00713
00714 }
00715
00716
00717
00718
00719
00720
00721
00722
00723
00724
00725
00726
00727
00728
00729
00730
00731
00732
00733
00734
00735
00736
00737
00738
00739 void
00740 spWhereSingular( eMatrix, pRow, pCol )
00741
00742 char *eMatrix;
00743 int *pRow, *pCol;
00744 {
00745 MatrixPtr Matrix = (MatrixPtr)eMatrix;
00746
00747
00748 ASSERT( IS_SPARSE( Matrix ) );
00749
00750 if (Matrix->Error == spSINGULAR OR Matrix->Error == spZERO_DIAG)
00751 { *pRow = Matrix->SingularRow;
00752 *pCol = Matrix->SingularCol;
00753 }
00754 else *pRow = *pCol = 0;
00755 return;
00756 }
00757
00758
00759
00760
00761
00762
00763
00764
00765
00766
00767
00768
00769
00770
00771
00772
00773
00774
00775
00776
00777
00778
00779 int
00780 spGetSize( eMatrix, External )
00781
00782 char *eMatrix;
00783 SPBOOLEAN External;
00784 {
00785 MatrixPtr Matrix = (MatrixPtr)eMatrix;
00786
00787
00788 ASSERT( IS_SPARSE( Matrix ) );
00789
00790 #if TRANSLATE
00791 if (External)
00792 return Matrix->ExtSize;
00793 else
00794 return Matrix->Size;
00795 #else
00796 return Matrix->Size;
00797 #endif
00798 }
00799
00800
00801
00802
00803
00804
00805
00806
00807
00808
00809
00810
00811
00812
00813
00814
00815
00816
00817 void
00818 spSetReal( eMatrix )
00819
00820 char *eMatrix;
00821 {
00822
00823
00824 ASSERT( IS_SPARSE( (MatrixPtr)eMatrix ) AND REAL);
00825 ((MatrixPtr)eMatrix)->Complex = NO;
00826 return;
00827 }
00828
00829
00830 void
00831 spSetComplex( eMatrix )
00832
00833 char *eMatrix;
00834 {
00835
00836
00837 ASSERT( IS_SPARSE( (MatrixPtr)eMatrix ) AND spCOMPLEX);
00838 ((MatrixPtr)eMatrix)->Complex = YES;
00839 return;
00840 }
00841
00842
00843
00844
00845
00846
00847
00848
00849
00850
00851
00852
00853
00854
00855
00856
00857
00858
00859
00860
00861 int
00862 spFillinCount( eMatrix )
00863
00864 char *eMatrix;
00865 {
00866
00867
00868 ASSERT( IS_SPARSE( (MatrixPtr)eMatrix ) );
00869 return ((MatrixPtr)eMatrix)->Fillins;
00870 }
00871
00872
00873 int
00874 spElementCount( eMatrix )
00875
00876 char *eMatrix;
00877 {
00878
00879
00880 ASSERT( IS_SPARSE( (MatrixPtr)eMatrix ) );
00881 return ((MatrixPtr)eMatrix)->Elements;
00882 }