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 #ifndef lint
00045
00046
00047 #endif
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061
00062
00063
00064
00065 #define spNO_MEMORY 4
00066
00067
00068
00069
00070
00071
00072
00073
00074
00075
00076
00077
00078
00079
00080 #define spINSIDE_SPARSE
00081 #include "spconfig.h"
00082 #include "spmatrix.h"
00083
00084 static void AllocateBlockOfAllocationList(MatrixPtr Matrix)
00085 {
00086 int I;
00087 AllocationListPtr ListPtr;
00088
00089
00090
00091 ListPtr = ALLOC(struct AllocationRecord, (ELEMENTS_PER_ALLOCATION+1));
00092 if (ListPtr == NULL)
00093 {
00094 Matrix->Error = spNO_MEMORY; return;
00095 }
00096
00097
00098
00099
00100 ListPtr->NextRecord = Matrix->TopOfAllocationList;
00101 Matrix->TopOfAllocationList = ListPtr;
00102 ListPtr += ELEMENTS_PER_ALLOCATION;
00103 for (I = ELEMENTS_PER_ALLOCATION; I > 0; I--)
00104 { ListPtr->NextRecord = ListPtr - 1;
00105 ListPtr--;
00106 }
00107
00108 Matrix->TopOfAllocationList->AllocatedPtr = (char *)ListPtr;
00109 Matrix->RecordsRemaining = ELEMENTS_PER_ALLOCATION;
00110 return;
00111 }
00112
00113 static void InitializeElementBlocks(MatrixPtr Matrix, int InitialNumberOfElements, int NumberOfFillinsExpected );
00114 static void RecordAllocation(MatrixPtr Matrix, char *AllocatedPtr);
00115
00116
00117
00118
00119
00120
00121
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
00151
00152 MatrixPtr spCreate(int Size, int *pError)
00153 {
00154 unsigned SizePlusOne;
00155 MatrixPtr Matrix;
00156 int I, AllocatedSize;
00157
00158 *pError = spOKAY;
00159
00160 AllocatedSize = MAX(Size, MINIMUM_ALLOCATED_SIZE);
00161 SizePlusOne = (unsigned)(AllocatedSize + 1);
00162
00163 if ((Matrix = ALLOC(struct MatrixFrame, 1)) == NULL)
00164 { *pError = spNO_MEMORY;
00165 return NULL;
00166 }
00167
00168
00169 Matrix->ID = SPARSE_ID;
00170 Matrix->Complex = 0;
00171 Matrix->PreviousMatrixWasComplex = 0;
00172 Matrix->Factored = 0;
00173 Matrix->Elements = 0;
00174 Matrix->Error = *pError;
00175 Matrix->Fillins = 0;
00176 Matrix->Reordered = 0;
00177 Matrix->NeedsOrdering = 1;
00178 Matrix->NumberOfInterchangesIsOdd = 0;
00179 Matrix->Partitioned = 0;
00180 Matrix->RowsLinked = 0;
00181 Matrix->InternalVectorsAllocated = 0;
00182 Matrix->SingularCol = 0;
00183 Matrix->SingularRow = 0;
00184 Matrix->Size = Size;
00185 Matrix->AllocatedSize = AllocatedSize;
00186 Matrix->ExtSize = Size;
00187 Matrix->AllocatedExtSize = AllocatedSize;
00188 Matrix->CurrentSize = 0;
00189 Matrix->ExtToIntColMap = NULL;
00190 Matrix->ExtToIntRowMap = NULL;
00191 Matrix->IntToExtColMap = NULL;
00192 Matrix->IntToExtRowMap = NULL;
00193 Matrix->MarkowitzRow = NULL;
00194 Matrix->MarkowitzCol = NULL;
00195 Matrix->MarkowitzProd = NULL;
00196 Matrix->DoCmplxDirect = NULL;
00197 Matrix->DoRealDirect = NULL;
00198 Matrix->Intermediate = NULL;
00199 Matrix->RelThreshold = DEFAULT_THRESHOLD;
00200 Matrix->AbsThreshold = 0.0;
00201
00202 Matrix->TopOfAllocationList = NULL;
00203 Matrix->RecordsRemaining = 0;
00204 Matrix->ElementsRemaining = 0;
00205 Matrix->FillinsRemaining = 0;
00206
00207 RecordAllocation( Matrix, (char *)Matrix );
00208 if (Matrix->Error == spNO_MEMORY) goto MemoryError;
00209
00210
00211 Matrix->TrashCan.Real = 0.0;
00212 Matrix->TrashCan.Row = 0;
00213 Matrix->TrashCan.Col = 0;
00214 Matrix->TrashCan.NextInRow = NULL;
00215 Matrix->TrashCan.NextInCol = NULL;
00216 #if INITIALIZE
00217 Matrix->TrashCan.pInitInfo = NULL;
00218 #endif
00219
00220
00221 CALLOC( Matrix->Diag, ElementPtr, SizePlusOne);
00222 if (Matrix->Diag == NULL)
00223 goto MemoryError;
00224
00225
00226 CALLOC( Matrix->FirstInCol, ElementPtr, SizePlusOne);
00227 if (Matrix->FirstInCol == NULL)
00228 goto MemoryError;
00229
00230
00231 CALLOC( Matrix->FirstInRow, ElementPtr, SizePlusOne);
00232 if (Matrix->FirstInRow == NULL)
00233 goto MemoryError;
00234
00235
00236 if (( Matrix->IntToExtColMap = ALLOC(int, SizePlusOne)) == NULL)
00237 goto MemoryError;
00238
00239
00240 if (( Matrix->IntToExtRowMap = ALLOC(int, SizePlusOne)) == NULL)
00241 goto MemoryError;
00242
00243
00244 for (I = 1; I <= AllocatedSize; I++)
00245 { Matrix->IntToExtRowMap[I] = I;
00246 Matrix->IntToExtColMap[I] = I;
00247 }
00248
00249 #if TRANSLATE
00250
00251 if (( Matrix->ExtToIntColMap = ALLOC(int, SizePlusOne)) == NULL)
00252 goto MemoryError;
00253
00254
00255 if (( Matrix->ExtToIntRowMap = ALLOC(int, SizePlusOne)) == NULL)
00256 goto MemoryError;
00257
00258
00259 for (I = 1; I <= AllocatedSize; I++)
00260 { Matrix->ExtToIntColMap[I] = -1;
00261 Matrix->ExtToIntRowMap[I] = -1;
00262 }
00263 Matrix->ExtToIntColMap[0] = 0;
00264 Matrix->ExtToIntRowMap[0] = 0;
00265 #endif
00266
00267
00268
00269 InitializeElementBlocks( Matrix, SPACE_FOR_ELEMENTS*AllocatedSize,
00270 SPACE_FOR_FILL_INS*AllocatedSize );
00271 if (Matrix->Error == spNO_MEMORY)
00272 goto MemoryError;
00273
00274 return Matrix;
00275
00276 MemoryError:
00277
00278
00279
00280 *pError = spNO_MEMORY;
00281 spDestroy(Matrix);
00282 return NULL;
00283 }
00284
00285
00286
00287
00288
00289
00290
00291
00292
00293
00294
00295
00296
00297
00298
00299
00300
00301
00302
00303
00304
00305
00306
00307 ElementPtr spcGetElement(MatrixPtr Matrix)
00308 {
00309 ElementPtr pElement;
00310
00311 if (Matrix->ElementsRemaining == 0)
00312 {
00313 pElement = ALLOC(struct MatrixElement, ELEMENTS_PER_ALLOCATION);
00314 RecordAllocation(Matrix, (char *)pElement);
00315 if (Matrix->Error == spNO_MEMORY) return NULL;
00316 Matrix->ElementsRemaining = ELEMENTS_PER_ALLOCATION;
00317 Matrix->NextAvailElement = pElement;
00318 }
00319
00320 Matrix->ElementsRemaining--;
00321 return Matrix->NextAvailElement++;
00322 }
00323
00324
00325
00326
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 static void InitializeElementBlocks(MatrixPtr Matrix, int InitialNumberOfElements, int NumberOfFillinsExpected )
00354 {
00355 ElementPtr pElement;
00356
00357
00358
00359
00360 pElement = ALLOC(struct MatrixElement, InitialNumberOfElements);
00361 RecordAllocation( Matrix, (char *)pElement );
00362 if (Matrix->Error == spNO_MEMORY) return;
00363 Matrix->ElementsRemaining = InitialNumberOfElements;
00364 Matrix->NextAvailElement = pElement;
00365
00366
00367 pElement = ALLOC(struct MatrixElement, NumberOfFillinsExpected);
00368 RecordAllocation( Matrix, (char *)pElement );
00369 if (Matrix->Error == spNO_MEMORY) return;
00370 Matrix->FillinsRemaining = NumberOfFillinsExpected;
00371 Matrix->NextAvailFillin = pElement;
00372
00373
00374 Matrix->FirstFillinListNode = ALLOC(struct FillinListNodeStruct,1);
00375 RecordAllocation( Matrix, (char *)Matrix->FirstFillinListNode );
00376 if (Matrix->Error == spNO_MEMORY) return;
00377 Matrix->LastFillinListNode = Matrix->FirstFillinListNode;
00378
00379 Matrix->FirstFillinListNode->pFillinList = pElement;
00380 Matrix->FirstFillinListNode->NumberOfFillinsInList =NumberOfFillinsExpected;
00381 Matrix->FirstFillinListNode->Next = NULL;
00382
00383 return;
00384 }
00385
00386
00387
00388
00389
00390
00391
00392
00393
00394
00395
00396
00397
00398
00399
00400
00401
00402
00403
00404 ElementPtr spcGetFillin(MatrixPtr Matrix)
00405 {
00406 struct FillinListNodeStruct *pListNode;
00407 ElementPtr pFillins;
00408
00409
00410
00411 #if !STRIP || LINT
00412 if (Matrix->FillinsRemaining == 0)
00413 return spcGetElement( Matrix );
00414 #endif
00415 #if STRIP || LINT
00416
00417 if (Matrix->FillinsRemaining == 0)
00418 { pListNode = Matrix->LastFillinListNode;
00419
00420
00421 if (pListNode->Next != NULL)
00422 { Matrix->LastFillinListNode = pListNode = pListNode->Next;
00423 Matrix->FillinsRemaining = pListNode->NumberOfFillinsInList;
00424 Matrix->NextAvailFillin = pListNode->pFillinList;
00425 }
00426 else
00427 {
00428
00429 pFillins = ALLOC(struct MatrixElement, ELEMENTS_PER_ALLOCATION);
00430 RecordAllocation( Matrix, (char *)pFillins );
00431 if (Matrix->Error == spNO_MEMORY) return NULL;
00432 Matrix->FillinsRemaining = ELEMENTS_PER_ALLOCATION;
00433 Matrix->NextAvailFillin = pFillins;
00434
00435
00436 pListNode->Next = ALLOC(struct FillinListNodeStruct,1);
00437 RecordAllocation( Matrix, (char *)pListNode->Next );
00438 if (Matrix->Error == spNO_MEMORY) return NULL;
00439 Matrix->LastFillinListNode = pListNode = pListNode->Next;
00440
00441 pListNode->pFillinList = pFillins;
00442 pListNode->NumberOfFillinsInList = ELEMENTS_PER_ALLOCATION;
00443 pListNode->Next = NULL;
00444 }
00445 }
00446 #endif
00447
00448
00449 Matrix->FillinsRemaining--;
00450 return Matrix->NextAvailFillin++;
00451 }
00452
00453
00454
00455
00456
00457
00458
00459
00460
00461
00462
00463
00464
00465
00466
00467
00468
00469 static void RecordAllocation(MatrixPtr Matrix, char *AllocatedPtr)
00470 {
00471
00472
00473
00474
00475 if (AllocatedPtr == NULL)
00476 {
00477 Matrix->Error = spNO_MEMORY; return;
00478 }
00479
00480 if (Matrix->RecordsRemaining == 0)
00481 {
00482 AllocateBlockOfAllocationList( Matrix );
00483 if (Matrix->Error == spNO_MEMORY)
00484 {
00485 FREE(AllocatedPtr);
00486 return;
00487 }
00488 }
00489
00490 (++Matrix->TopOfAllocationList)->AllocatedPtr = AllocatedPtr;
00491 Matrix->RecordsRemaining--;
00492 return;
00493 }
00494
00495
00496
00497
00498
00499
00500
00501
00502
00503
00504
00505
00506
00507
00508
00509
00510
00511
00512
00513
00514 void spDestroy(MatrixPtr Matrix)
00515 {
00516 AllocationListPtr ListPtr, NextListPtr;
00517
00518
00519 ASSERT( IS_SPARSE( Matrix ) );
00520
00521
00522 FREE( Matrix->IntToExtColMap );
00523 FREE( Matrix->IntToExtRowMap );
00524 FREE( Matrix->ExtToIntColMap );
00525 FREE( Matrix->ExtToIntRowMap );
00526 FREE( Matrix->Diag );
00527 FREE( Matrix->FirstInRow );
00528 FREE( Matrix->FirstInCol );
00529 FREE( Matrix->MarkowitzRow );
00530 FREE( Matrix->MarkowitzCol );
00531 FREE( Matrix->MarkowitzProd );
00532 FREE( Matrix->DoCmplxDirect );
00533 FREE( Matrix->DoRealDirect );
00534 FREE( Matrix->Intermediate );
00535
00536
00537 ListPtr = Matrix->TopOfAllocationList;
00538 while (ListPtr != NULL)
00539 {
00540 NextListPtr = ListPtr->NextRecord;
00541 FREE( ListPtr->AllocatedPtr );
00542 ListPtr = NextListPtr;
00543 }
00544 return;
00545 }
00546
00547
00548
00549
00550
00551
00552
00553
00554
00555
00556
00557
00558 int spError(MatrixPtr Matrix )
00559 {
00560 if (Matrix != NULL)
00561 {
00562 ASSERT(Matrix->ID == SPARSE_ID); return Matrix->Error;
00563 }
00564 else return spNO_MEMORY;
00565 }
00566
00567
00568
00569
00570
00571
00572
00573
00574
00575
00576
00577
00578
00579
00580
00581 void spWhereSingular(MatrixPtr Matrix, int* pRow, int* pCol )
00582 {
00583 ASSERT( IS_SPARSE( Matrix ) );
00584
00585 if (Matrix->Error == spSINGULAR || Matrix->Error == spZERO_DIAG)
00586 { *pRow = Matrix->SingularRow;
00587 *pCol = Matrix->SingularCol;
00588 }
00589 else *pRow = *pCol = 0;
00590 return;
00591 }
00592
00593
00594
00595
00596
00597
00598
00599
00600
00601
00602
00603
00604
00605
00606
00607 int spGetSize(MatrixPtr Matrix, BOOLEAN External)
00608 {
00609 ASSERT( IS_SPARSE( Matrix ) );
00610 #if TRANSLATE
00611 if (External)
00612 return Matrix->ExtSize;
00613 else
00614 return Matrix->Size;
00615 #else
00616 return Matrix->Size;
00617 #endif
00618 }
00619
00620 int spElementCount(MatrixPtr Matrix)
00621 {
00622 ASSERT( IS_SPARSE( Matrix ) );
00623 return Matrix->Elements;
00624 }