00001 /* 00002 * MATRIX SOLVE MODULE 00003 * 00004 * Author: Advising professor: 00005 * Kenneth S. Kundert Alberto Sangiovanni-Vincentelli 00006 * UC Berkeley 00007 * 00008 * This file contains the forward and backward substitution routines for 00009 * the sparse matrix routines. 00010 * 00011 * >>> User accessible functions contained in this file: 00012 * spSolve 00013 * spSolveTransposed 00014 * 00015 * >>> Other functions contained in this file: 00016 * SolveComplexMatrix 00017 * SolveComplexTransposedMatrix 00018 */ 00019 00020 00021 /* 00022 * Revision and copyright information. 00023 * 00024 * Copyright (c) 1985,86,87,88 00025 * by Kenneth S. Kundert and the University of California. 00026 * 00027 * Permission to use, copy, modify, and distribute this software and 00028 * its documentation for any purpose and without fee is hereby granted, 00029 * provided that the copyright notices appear in all copies and 00030 * supporting documentation and that the authors and the University of 00031 * California are properly credited. The authors and the University of 00032 * California make no representations as to the suitability of this 00033 * software for any purpose. It is provided `as is', without express 00034 * or implied warranty. 00035 */ 00036 00037 /* 00038 * IMPORTS 00039 * 00040 * >>> Import descriptions: 00041 * spConfig.h 00042 * Macros that customize the sparse matrix routines. 00043 * spmatrix.h 00044 * Macros and declarations to be imported by the user. 00045 * spDefs.h 00046 * Matrix type and macro definitions for the sparse matrix routines. 00047 */ 00048 00049 #define spINSIDE_SPARSE 00050 #include "spconfig.h" 00051 #include "spmatrix.h" 00052 00053 #if spCOMPLEX 00054 /* 00055 * SOLVE COMPLEX MATRIX EQUATION 00056 * 00057 * Performs forward elimination and back substitution to find the 00058 * unknown vector from the RHS vector and factored matrix. This 00059 * routine assumes that the pivots are associated with the lower 00060 * triangular (L) matrix and that the diagonal of the upper triangular 00061 * (U) matrix consists of ones. This routine arranges the computation 00062 * in different way than is traditionally used in order to exploit the 00063 * sparsity of the right-hand side. See the reference in spRevision. 00064 * 00065 * >>> Arguments: 00066 * Matrix <input> (char *) 00067 * Pointer to matrix. 00068 * RHS <input> (double*) 00069 * RHS is the real portion of the input data array, the right hand 00070 * side. This data is undisturbed and may be reused for other solves. 00071 * Solution <output> (double*) 00072 * Solution is the real portion of the output data array. This routine 00073 * is constructed such that RHS and Solution can be the same 00074 * array. 00075 * iRHS <input> (double*) 00076 * iRHS is the imaginary portion of the input data array, the right 00077 * hand side. This data is undisturbed and may be reused for other solves. 00078 * If spSEPARATED_COMPLEX_VECTOR is set false, there is no need to 00079 * supply this array. 00080 * iSolution <output> (double*) 00081 * iSolution is the imaginary portion of the output data array. This 00082 * routine is constructed such that iRHS and iSolution can be 00083 * the same array. If spSEPARATED_COMPLEX_VECTOR is set false, there is no 00084 * need to supply this array. 00085 * 00086 * >>> Local variables: 00087 * Intermediate (ComplexVector) 00088 * Temporary storage for use in forward elimination and backward 00089 * substitution. Commonly referred to as c, when the LU factorization 00090 * equations are given as Ax = b, Lc = b, Ux = c. 00091 * Local version of Matrix->Intermediate, which was created during 00092 * the initial factorization in function CreateInternalVectors() in the 00093 * matrix factorization module. 00094 * pElement (ElementPtr) 00095 * Pointer used to address elements in both the lower and upper triangle 00096 * matrices. 00097 * pExtOrder (int *) 00098 * Pointer used to sequentially access each entry in IntToExtRowMap 00099 * and IntToExtColMap arrays. Used to quickly scramble and unscramble 00100 * RHS and Solution to account for row and column interchanges. 00101 * pPivot (ElementPtr) 00102 * Pointer that points to current pivot or diagonal element. 00103 * Size (int) 00104 * Size of matrix. Made local to reduce indirection. 00105 * Temp (ComplexNumber) 00106 * Temporary storage for entries in arrays. 00107 * 00108 * >>> Obscure Macros 00109 * IMAG_VECTORS 00110 * Replaces itself with `, iRHS, iSolution' if the options spCOMPLEX and 00111 * spSEPARATED_COMPLEX_VECTORS are set, otherwise it disappears 00112 * without a trace. 00113 */ 00114 00115 static void SolveComplexMatrix(MatrixPtr Matrix, double *RHS, double *Solution IMAG_VECTORS) 00116 { 00117 ElementPtr pElement; 00118 ComplexVector Intermediate; 00119 int I, *pExtOrder, Size; 00120 ElementPtr pPivot; 00121 ComplexVector ExtVector; 00122 ComplexNumber Temp; 00123 Size = Matrix->Size; 00124 Intermediate = (ComplexVector)Matrix->Intermediate; 00125 00126 /* Correct array pointers for ARRAY_OFFSET. */ 00127 #if !ARRAY_OFFSET 00128 #if spSEPARATED_COMPLEX_VECTORS 00129 --RHS; --iRHS; 00130 --Solution; --iSolution; 00131 #else 00132 RHS -= 2; Solution -= 2; 00133 #endif 00134 #endif 00135 00136 /* Initialize Intermediate vector. */ 00137 pExtOrder = &Matrix->IntToExtRowMap[Size]; 00138 00139 #if spSEPARATED_COMPLEX_VECTORS 00140 for (I = Size; I > 0; I--) 00141 { Intermediate[I].Real = RHS[*(pExtOrder)]; 00142 Intermediate[I].Imag = iRHS[*(pExtOrder--)]; 00143 } 00144 #else 00145 ExtVector = (ComplexVector)RHS; 00146 for (I = Size; I > 0; I--) 00147 Intermediate[I] = ExtVector[*(pExtOrder--)]; 00148 #endif 00149 00150 /* Forward substitution. Solves Lc = b.*/ 00151 for (I = 1; I <= Size; I++) 00152 { 00153 Temp = Intermediate[I]; 00154 00155 /* This step of the substitution is skipped if Temp equals zero. */ 00156 if ((Temp.Real != 0.0) || (Temp.Imag != 0.0)) 00157 { 00158 pPivot = Matrix->Diag[I]; 00159 /* Cmplx expr: Temp *= (1.0 / Pivot). */ 00160 CMPLX_MULT_ASSIGN(Temp, *pPivot); 00161 Intermediate[I] = Temp; 00162 pElement = pPivot->NextInCol; 00163 while (pElement != NULL) 00164 { 00165 /* Cmplx expr: Intermediate[Element->Row] -= Temp * *Element. */ 00166 CMPLX_MULT_SUBT_ASSIGN(Intermediate[pElement->Row], 00167 Temp, *pElement); 00168 pElement = pElement->NextInCol; 00169 } 00170 } 00171 } 00172 00173 /* Backward Substitution. Solves Ux = c.*/ 00174 for (I = Size; I > 0; I--) 00175 { 00176 Temp = Intermediate[I]; 00177 pElement = Matrix->Diag[I]->NextInRow; 00178 00179 while (pElement != NULL) 00180 { 00181 /* Cmplx expr: Temp -= *Element * Intermediate[Element->Col]. */ 00182 CMPLX_MULT_SUBT_ASSIGN(Temp, *pElement,Intermediate[pElement->Col]); 00183 pElement = pElement->NextInRow; 00184 } 00185 Intermediate[I] = Temp; 00186 } 00187 00188 /* Unscramble Intermediate vector while placing data in to Solution vector. */ 00189 pExtOrder = &Matrix->IntToExtColMap[Size]; 00190 00191 #if spSEPARATED_COMPLEX_VECTORS 00192 for (I = Size; I > 0; I--) 00193 { 00194 Solution[*(pExtOrder)] = Intermediate[I].Real; 00195 iSolution[*(pExtOrder--)] = Intermediate[I].Imag; 00196 } 00197 #else 00198 ExtVector = (ComplexVector)Solution; 00199 for (I = Size; I > 0; I--) 00200 ExtVector[*(pExtOrder--)] = Intermediate[I]; 00201 #endif 00202 00203 return; 00204 } 00205 #endif /* spCOMPLEX */ 00206 #if TRANSPOSE && spCOMPLEX 00207 /* 00208 * SOLVE COMPLEX TRANSPOSED MATRIX EQUATION 00209 * 00210 * Performs forward elimination and back substitution to find the 00211 * unknown vector from the RHS vector and transposed factored 00212 * matrix. This routine is useful when performing sensitivity analysis 00213 * on a circuit using the adjoint method. This routine assumes that 00214 * the pivots are associated with the untransposed lower triangular 00215 * (L) matrix and that the diagonal of the untransposed upper 00216 * triangular (U) matrix consists of ones. 00217 * 00218 * >>> Arguments: 00219 * Matrix <input> (char *) 00220 * Pointer to matrix. 00221 * RHS <input> (double*) 00222 * RHS is the input data array, the right hand 00223 * side. This data is undisturbed and may be reused for other solves. 00224 * This vector is only the real portion if the matrix is complex and 00225 * spSEPARATED_COMPLEX_VECTORS is set true. 00226 * Solution <output> (double*) 00227 * Solution is the real portion of the output data array. This routine 00228 * is constructed such that RHS and Solution can be the same array. 00229 * This vector is only the real portion if the matrix is complex and 00230 * spSEPARATED_COMPLEX_VECTORS is set true. 00231 * iRHS <input> (double*) 00232 * iRHS is the imaginary portion of the input data array, the right 00233 * hand side. This data is undisturbed and may be reused for other solves. 00234 * If either spCOMPLEX or spSEPARATED_COMPLEX_VECTOR is set false, there 00235 * is no need to supply this array. 00236 * iSolution <output> (double*) 00237 * iSolution is the imaginary portion of the output data array. This 00238 * routine is constructed such that iRHS and iSolution can be 00239 * the same array. If spCOMPLEX or spSEPARATED_COMPLEX_VECTOR is set 00240 * false, there is no need to supply this array. 00241 * 00242 * >>> Local variables: 00243 * Intermediate (ComplexVector) 00244 * Temporary storage for use in forward elimination and backward 00245 * substitution. Commonly referred to as c, when the LU factorization 00246 * equations are given as Ax = b, Lc = b, Ux = c. Local version of 00247 * Matrix->Intermediate, which was created during 00248 * the initial factorization in function CreateInternalVectors() in the 00249 * matrix factorization module. 00250 * pElement (ElementPtr) 00251 * Pointer used to address elements in both the lower and upper triangle 00252 * matrices. 00253 * pExtOrder (int *) 00254 * Pointer used to sequentially access each entry in IntToExtRowMap 00255 * and IntToExtColMap arrays. Used to quickly scramble and unscramble 00256 * RHS and Solution to account for row and column interchanges. 00257 * pPivot (ElementPtr) 00258 * Pointer that points to current pivot or diagonal element. 00259 * Size (int) 00260 * Size of matrix. Made local to reduce indirection. 00261 * Temp (ComplexNumber) 00262 * Temporary storage for entries in arrays. 00263 * 00264 * >>> Obscure Macros 00265 * IMAG_VECTORS 00266 * Replaces itself with `, iRHS, iSolution' if the options spCOMPLEX and 00267 * spSEPARATED_COMPLEX_VECTORS are set, otherwise it disappears 00268 * without a trace. 00269 */ 00270 00271 static void SolveComplexTransposedMatrix(MatrixPtr Matrix, double *RHS, double *Solution IMAG_VECTORS) 00272 { 00273 ElementPtr pElement; 00274 ComplexVector Intermediate; 00275 int I, *pExtOrder, Size; 00276 ComplexVector ExtVector; 00277 ElementPtr pPivot; 00278 ComplexNumber Temp; 00279 Size = Matrix->Size; 00280 Intermediate = (ComplexVector)Matrix->Intermediate; 00281 00282 /* Correct array pointers for ARRAY_OFFSET. */ 00283 #if !ARRAY_OFFSET 00284 #if spSEPARATED_COMPLEX_VECTORS 00285 --RHS; --iRHS; --Solution; --iSolution; 00286 #else 00287 RHS -= 2; Solution -= 2; 00288 #endif 00289 #endif 00290 00291 /* Initialize Intermediate vector. */ 00292 pExtOrder = &Matrix->IntToExtColMap[Size]; 00293 00294 #if spSEPARATED_COMPLEX_VECTORS 00295 for (I = Size; I > 0; I--) 00296 { 00297 Intermediate[I].Real = RHS[*(pExtOrder)]; 00298 Intermediate[I].Imag = iRHS[*(pExtOrder--)]; 00299 } 00300 #else 00301 ExtVector = (ComplexVector)RHS; 00302 for (I = Size; I > 0; I--) 00303 Intermediate[I] = ExtVector[*(pExtOrder--)]; 00304 #endif 00305 00306 /* Forward elimination. */ 00307 for (I = 1; I <= Size; I++) 00308 { 00309 Temp = Intermediate[I]; 00310 00311 /* This step of the elimination is skipped if Temp equals zero. */ 00312 if ((Temp.Real != 0.0) || (Temp.Imag != 0.0)) 00313 { 00314 pElement = Matrix->Diag[I]->NextInRow; 00315 while (pElement != NULL) 00316 { 00317 /* Cmplx expr: Intermediate[Element->Col] -= Temp * *Element. */ 00318 CMPLX_MULT_SUBT_ASSIGN( Intermediate[pElement->Col], Temp, *pElement); 00319 pElement = pElement->NextInRow; 00320 } 00321 } 00322 } 00323 00324 /* Backward Substitution. */ 00325 for (I = Size; I > 0; I--) 00326 { 00327 pPivot = Matrix->Diag[I]; 00328 Temp = Intermediate[I]; 00329 pElement = pPivot->NextInCol; 00330 while (pElement != NULL) 00331 { 00332 /* Cmplx expr: Temp -= Intermediate[Element->Row] * *Element. */ 00333 CMPLX_MULT_SUBT_ASSIGN(Temp,Intermediate[pElement->Row],*pElement); 00334 pElement = pElement->NextInCol; 00335 } 00336 /* Cmplx expr: Intermediate = Temp * (1.0 / *pPivot). */ 00337 CMPLX_MULT(Intermediate[I], Temp, *pPivot); 00338 } 00339 00340 /* Unscramble Intermediate vector while placing data in to Solution vector. */ 00341 pExtOrder = &Matrix->IntToExtRowMap[Size]; 00342 00343 #if spSEPARATED_COMPLEX_VECTORS 00344 for (I = Size; I > 0; I--) 00345 { 00346 Solution[*(pExtOrder)] = Intermediate[I].Real; 00347 iSolution[*(pExtOrder--)] = Intermediate[I].Imag; 00348 } 00349 #else 00350 ExtVector = (ComplexVector)Solution; 00351 for (I = Size; I > 0; I--) 00352 ExtVector[*(pExtOrder--)] = Intermediate[I]; 00353 #endif 00354 return; 00355 } 00356 #endif /* TRANSPOSE && spCOMPLEX */ 00357 00358 /* 00359 * SOLVE MATRIX EQUATION 00360 * 00361 * Performs forward elimination and back substitution to find the 00362 * unknown vector from the RHS vector and factored matrix. This 00363 * routine assumes that the pivots are associated with the lower 00364 * triangular (L) matrix and that the diagonal of the upper triangular 00365 * (U) matrix consists of ones. This routine arranges the computation 00366 * in different way than is traditionally used in order to exploit the 00367 * sparsity of the right-hand side. See the reference in spRevision. 00368 * 00369 * >>> Arguments: 00370 * Matrix <input> (char *) 00371 * Pointer to matrix. 00372 * RHS <input> (double*) 00373 * RHS is the input data array, the right hand side. This data is 00374 * undisturbed and may be reused for other solves. 00375 * Solution <output> (double*) 00376 * Solution is the output data array. This routine is constructed such that 00377 * RHS and Solution can be the same array. 00378 * iRHS <input> (double*) 00379 * iRHS is the imaginary portion of the input data array, the right 00380 * hand side. This data is undisturbed and may be reused for other solves. 00381 * This argument is only necessary if matrix is complex and if 00382 * spSEPARATED_COMPLEX_VECTOR is set true. 00383 * iSolution <output> (double*) 00384 * iSolution is the imaginary portion of the output data array. This 00385 * routine is constructed such that iRHS and iSolution can be 00386 * the same array. This argument is only necessary if matrix is complex 00387 * and if spSEPARATED_COMPLEX_VECTOR is set true. 00388 * 00389 * >>> Local variables: 00390 * Intermediate (double*) 00391 * Temporary storage for use in forward elimination and backward 00392 * substitution. Commonly referred to as c, when the LU factorization 00393 * equations are given as Ax = b, Lc = b, Ux = c Local version of 00394 * Matrix->Intermediate, which was created during the initial 00395 * factorization in function CreateInternalVectors() in the matrix 00396 * factorization module. 00397 * pElement (ElementPtr) 00398 * Pointer used to address elements in both the lower and upper triangle 00399 * matrices. 00400 * pExtOrder (int *) 00401 * Pointer used to sequentially access each entry in IntToExtRowMap 00402 * and IntToExtColMap arrays. Used to quickly scramble and unscramble 00403 * RHS and Solution to account for row and column interchanges. 00404 * pPivot (ElementPtr) 00405 * Pointer that points to current pivot or diagonal element. 00406 * Size (int) 00407 * Size of matrix. Made local to reduce indirection. 00408 * Temp (double) 00409 * Temporary storage for entries in arrays. 00410 * 00411 * >>> Obscure Macros 00412 * IMAG_VECTORS 00413 * Replaces itself with `, iRHS, iSolution' if the options spCOMPLEX and 00414 * spSEPARATED_COMPLEX_VECTORS are set, otherwise it disappears 00415 * without a trace. 00416 */ 00417 00418 void spSolve(MatrixPtr Matrix, double *RHS, double *Solution) 00419 { 00420 ElementPtr pElement, pPivot; 00421 double* Intermediate, Temp; 00422 int I, *pExtOrder, Size; 00423 00424 //ASSERT( IS_VALID(Matrix) && IS_FACTORED(Matrix) ); 00425 00426 #if spCOMPLEX 00427 if (Matrix->Complex) 00428 { 00429 SolveComplexMatrix( Matrix, RHS, Solution IMAG_VECTORS ); 00430 return; 00431 } 00432 #endif 00433 00434 #if REAL 00435 Intermediate = Matrix->Intermediate; 00436 Size = Matrix->Size; 00437 00438 /* Correct array pointers for ARRAY_OFFSET. */ 00439 #if !ARRAY_OFFSET 00440 --RHS; --Solution; 00441 #endif 00442 00443 /* Initialize Intermediate vector. */ 00444 pExtOrder = &Matrix->IntToExtRowMap[Size]; 00445 for (I = Size; I > 0; I--) 00446 Intermediate[I] = RHS[*(pExtOrder--)]; 00447 00448 /* Forward elimination. Solves Lc = b.*/ 00449 for (I = 1; I <= Size; I++) 00450 { 00451 /* This step of the elimination is skipped if Temp equals zero. */ 00452 if ((Temp = Intermediate[I]) != 0.0) 00453 { 00454 pPivot = Matrix->Diag[I]; 00455 Intermediate[I] = (Temp *= pPivot->Real); 00456 pElement = pPivot->NextInCol; 00457 while (pElement != NULL) 00458 { 00459 Intermediate[pElement->Row] -= Temp * pElement->Real; 00460 pElement = pElement->NextInCol; 00461 } 00462 } 00463 } 00464 /* Backward Substitution. Solves Ux = c.*/ 00465 for (I = Size; I > 0; I--) 00466 { 00467 Temp = Intermediate[I]; 00468 pElement = Matrix->Diag[I]->NextInRow; 00469 while (pElement != NULL) 00470 { 00471 Temp -= pElement->Real * Intermediate[pElement->Col]; 00472 pElement = pElement->NextInRow; 00473 } 00474 Intermediate[I] = Temp; 00475 } 00476 00477 /* Unscramble Intermediate vector while placing data in to Solution vector. */ 00478 pExtOrder = &Matrix->IntToExtColMap[Size]; 00479 for (I = Size; I > 0; I--) 00480 Solution[*(pExtOrder--)] = Intermediate[I]; 00481 00482 return; 00483 #endif /* REAL */ 00484 } 00485 #if TRANSPOSE 00486 /* 00487 * SOLVE TRANSPOSED MATRIX EQUATION 00488 * 00489 * Performs forward elimination and back substitution to find the 00490 * unknown vector from the RHS vector and transposed factored 00491 * matrix. This routine is useful when performing sensitivity analysis 00492 * on a circuit using the adjoint method. This routine assumes that 00493 * the pivots are associated with the untransposed lower triangular 00494 * (L) matrix and that the diagonal of the untransposed upper 00495 * triangular (U) matrix consists of ones. 00496 * 00497 * >>> Arguments: 00498 * Matrix <input> (char *) 00499 * Pointer to matrix. 00500 * RHS <input> (double*) 00501 * RHS is the input data array, the right hand side. This data is 00502 * undisturbed and may be reused for other solves. 00503 * Solution <output> (double*) 00504 * Solution is the output data array. This routine is constructed such that 00505 * RHS and Solution can be the same array. 00506 * iRHS <input> (double*) 00507 * iRHS is the imaginary portion of the input data array, the right 00508 * hand side. This data is undisturbed and may be reused for other solves. 00509 * If spSEPARATED_COMPLEX_VECTOR is set false, or if matrix is real, there 00510 * is no need to supply this array. 00511 * iSolution <output> (double*) 00512 * iSolution is the imaginary portion of the output data array. This 00513 * routine is constructed such that iRHS and iSolution can be 00514 * the same array. If spSEPARATED_COMPLEX_VECTOR is set false, or if 00515 * matrix is real, there is no need to supply this array. 00516 * 00517 * >>> Local variables: 00518 * Intermediate (double*) 00519 * Temporary storage for use in forward elimination and backward 00520 * substitution. Commonly referred to as c, when the LU factorization 00521 * equations are given as Ax = b, Lc = b, Ux = c. Local version of 00522 * Matrix->Intermediate, which was created during the initial 00523 * factorization in function CreateInternalVectors() in the matrix 00524 * factorization module. 00525 * pElement (ElementPtr) 00526 * Pointer used to address elements in both the lower and upper triangle 00527 * matrices. 00528 * pExtOrder (int *) 00529 * Pointer used to sequentially access each entry in IntToExtRowMap 00530 * and IntToExtRowMap arrays. Used to quickly scramble and unscramble 00531 * RHS and Solution to account for row and column interchanges. 00532 * pPivot (ElementPtr) 00533 * Pointer that points to current pivot or diagonal element. 00534 * Size (int) 00535 * Size of matrix. Made local to reduce indirection. 00536 * Temp (double) 00537 * Temporary storage for entries in arrays. 00538 * 00539 * >>> Obscure Macros 00540 * IMAG_VECTORS 00541 * Replaces itself with `, iRHS, iSolution' if the options spCOMPLEX and 00542 * spSEPARATED_COMPLEX_VECTORS are set, otherwise it disappears 00543 * without a trace. 00544 */ 00545 00546 /*VARARGS3*/ 00547 00548 void spSolveTransposed(MatrixPtr Matrix, double *RHS, double *Solution) 00549 { 00550 ElementPtr pElement, pPivot; 00551 double* Intermediate, Temp; 00552 int I, *pExtOrder, Size; 00553 ASSERT( IS_VALID(Matrix) && IS_FACTORED(Matrix) ); 00554 00555 #if spCOMPLEX 00556 if (Matrix->Complex) 00557 { 00558 SolveComplexTransposedMatrix( Matrix, RHS, Solution IMAG_VECTORS ); 00559 return; 00560 } 00561 #endif 00562 00563 #if REAL 00564 Size = Matrix->Size; 00565 Intermediate = Matrix->Intermediate; 00566 00567 /* Correct array pointers for ARRAY_OFFSET. */ 00568 #if !ARRAY_OFFSET 00569 --RHS; --Solution; 00570 #endif 00571 00572 /* Initialize Intermediate vector. */ 00573 pExtOrder = &Matrix->IntToExtColMap[Size]; 00574 for (I = Size; I > 0; I--) 00575 Intermediate[I] = RHS[*(pExtOrder--)]; 00576 00577 /* Forward elimination. */ 00578 for (I = 1; I <= Size; I++) 00579 { 00580 /* This step of the elimination is skipped if Temp equals zero. */ 00581 if ((Temp = Intermediate[I]) != 0.0) 00582 { 00583 pElement = Matrix->Diag[I]->NextInRow; 00584 while (pElement != NULL) 00585 { 00586 Intermediate[pElement->Col] -= Temp * pElement->Real; 00587 pElement = pElement->NextInRow; 00588 } 00589 00590 } 00591 } 00592 00593 /* Backward Substitution. */ 00594 for (I = Size; I > 0; I--) 00595 { pPivot = Matrix->Diag[I]; 00596 Temp = Intermediate[I]; 00597 pElement = pPivot->NextInCol; 00598 while (pElement != NULL) 00599 { 00600 Temp -= pElement->Real * Intermediate[pElement->Row]; 00601 pElement = pElement->NextInCol; 00602 } 00603 Intermediate[I] = Temp * pPivot->Real; 00604 } 00605 00606 /* Unscramble Intermediate vector while placing data in to Solution vector. */ 00607 pExtOrder = &Matrix->IntToExtRowMap[Size]; 00608 for (I = Size; I > 0; I--) 00609 Solution[*(pExtOrder--)] = Intermediate[I]; 00610 return; 00611 #endif /* REAL */ 00612 } 00613 #endif /* TRANSPOSE */