00001 /* 00002 * SPARSE FORTRAN MODULE 00003 * 00004 * Author: Advising professor: 00005 * Kenneth S. Kundert Alberto Sangiovanni-Vincentelli 00006 * UC Berkeley 00007 * 00008 * This module contains routines that interface Sparse1.3 to a calling 00009 * program written in fortran. Almost every externally available Sparse1.3 00010 * routine has a counterpart defined in this file, with the name the 00011 * same except the `sp' prefix is changed to `sf'. The spADD_ELEMENT 00012 * and spADD_QUAD macros are also replaced with the sfAdd1 and sfAdd4 00013 * functions defined in this file. 00014 * 00015 * To ease porting this file to different operating systems, the names of 00016 * the functions can be easily redefined (search for `Routine Renaming'). 00017 * A simple example of a FORTRAN program that calls Sparse is included in 00018 * this file (search for Example). When interfacing to a FORTRAN program, 00019 * the ARRAY_OFFSET option should be set to NO (see spConfig.h). 00020 * 00021 * DISCLAIMER: 00022 * These interface routines were written by a C programmer who has little 00023 * experience with FORTRAN. The routines have had minimal testing. 00024 * Any interface between two languages is going to have portability 00025 * problems, this one is no exception. 00026 * 00027 * 00028 * >>> User accessible functions contained in this file: 00029 * sfCreate() 00030 * sfDestroy() 00031 * sfStripFills() 00032 * sfClear() 00033 * sfGetElement() 00034 * sfGetAdmittance() 00035 * sfGetQuad() 00036 * sfGetOnes() 00037 * sfAdd1Real() 00038 * sfAdd1Imag() 00039 * sfAdd1Complex() 00040 * sfAdd4Real() 00041 * sfAdd4Imag() 00042 * sfAdd4Complex() 00043 * sfOrderAndFactor() 00044 * sfFactor() 00045 * sfPartition() 00046 * sfSolve() 00047 * sfSolveTransposed() 00048 * sfPrint() 00049 * sfFileMatrix() 00050 * sfFileVector() 00051 * sfFileStats() 00052 * sfMNA_Preorder() 00053 * sfScale() 00054 * sfMultiply() 00055 * sfTransMultiply() 00056 * sfDeterminant() 00057 * sfError() 00058 * sfWhereSingular() 00059 * sfGetSize() 00060 * sfSetReal() 00061 * sfSetComplex() 00062 * sfFillinCount() 00063 * sfElementCount() 00064 * sfDeleteRowAndCol() 00065 * sfPseudoCondition() 00066 * sfCondition() 00067 * sfNorm() 00068 * sfLargestElement() 00069 * sfRoundoff() 00070 */ 00071 00072 /* 00073 * FORTRAN -- C COMPATIBILITY 00074 * 00075 * Fortran and C data types correspond in the following way: 00076 * -- C -- -- FORTRAN -- 00077 * int INTEGER*4 or INTEGER*2 (machine dependent, usually int*4) 00078 * long INTEGER*4 00079 * float REAL 00080 * double DOUBLE PRECISION (used by default in preference to float) 00081 * 00082 * The complex number format used by Sparse is compatible with that 00083 * used by FORTRAN. C pointers are passed to FORTRAN as longs, they should 00084 * not be used in any way in FORTRAN. 00085 */ 00086 00087 00088 00089 /* 00090 * Revision and copyright information. 00091 * 00092 * Copyright (c) 1985,86,87,88 00093 * by Kenneth S. Kundert and the University of California. 00094 * 00095 * Permission to use, copy, modify, and distribute this software and its 00096 * documentation for any purpose and without fee is hereby granted, provided 00097 * that the above copyright notice appear in all copies and supporting 00098 * documentation and that the authors and the University of California 00099 * are properly credited. The authors and the University of California 00100 * make no representations as to the suitability of this software for 00101 * any purpose. It is provided `as is', without express or implied warranty. 00102 */ 00103 00104 #ifndef lint 00105 static char copyright[] = 00106 "Sparse1.3: Copyright (c) 1985,86,87,88 by Kenneth S. Kundert"; 00107 static char RCSid[] = 00108 "@(#)$Header: spFortran.c,v 1.1 88/06/18 11:15:30 kundert Exp $"; 00109 #endif 00110 00111 00112 00113 00114 /* 00115 * IMPORTS 00116 * 00117 * >>> Import descriptions: 00118 * spConfig.h 00119 * Macros that customize the sparse matrix routines. 00120 * spmatrix.h 00121 * Macros and declarations to be imported by the user. 00122 * spDefs.h 00123 * Matrix type and macro definitions for the sparse matrix routines. 00124 */ 00125 00126 #define spINSIDE_SPARSE 00127 #include "spConfig.h" 00128 #include "spmatrix.h" 00129 #include "spDefs.h" 00130 00131 #if FORTRAN 00132 00133 00134 00135 00136 /* 00137 * Routine Renaming 00138 */ 00139 00140 #define sfCreate sfcreate 00141 #define sfStripFills sfstripfills 00142 #define sfDestroy sfdestroy 00143 #define sfClear sfzero 00144 #define sfGetElement sfgetelement 00145 #define sfGetAdmittance sfgetadmittance 00146 #define sfGetQuad sfgetquad 00147 #define sfGetOnes sfgetones 00148 #define sfAdd1Real sfadd1real 00149 #define sfAdd1Imag sfadd1imag 00150 #define sfAdd1Complex sfadd1complex 00151 #define sfAdd4Real sfadd4real 00152 #define sfAdd4Imag sfadd4imag 00153 #define sfAdd4Complex sfadd4complex 00154 #define sfOrderAndFactor sforderandfactor 00155 #define sfFactor sffactor 00156 #define sfPartition sfpartition 00157 #define sfSolve sfsolve 00158 #define sfSolveTransposed sfsolvetransposed 00159 #define sfPrint sfprint 00160 #define sfFileMatrix sffilematrix 00161 #define sfFileVector sffilevector 00162 #define sfFileStats sffilestats 00163 #define sfMNA_Preorder sfmna_preorder 00164 #define sfScale sfscale 00165 #define sfMultiply sfmultiply 00166 #define sfDeterminant sfdeterminant 00167 #define sfError sferror 00168 #define sfWhereSingular sfwheresingular 00169 #define sfGetSize sfgetsize 00170 #define sfSetReal sfsetreal 00171 #define sfSetComplex sfsetcomplex 00172 #define sfFillinCount sffillincount 00173 #define sfElementCount sfelementcount 00174 #define sfDeleteRowAndCol sfdeleterowandcol 00175 #define sfPseudoCondition sfpseudocondition 00176 #define sfCondition sfcondition 00177 #define sfNorm sfnorm 00178 #define sfLargestElement sflargestelement 00179 #define sfRoundoff sfroundoff 00180 00181 00182 00183 /* 00184 * Example of a FORTRAN Program Calling Sparse 00185 * 00186 00187 integer matrix, error, sfCreate, sfGetElement, spFactor 00188 integer element(10) 00189 double precision rhs(4), solution(4) 00190 c 00191 matrix = sfCreate(4,0,error) 00192 element(1) = sfGetElement(matrix,1,1) 00193 element(2) = sfGetElement(matrix,1,2) 00194 element(3) = sfGetElement(matrix,2,1) 00195 element(4) = sfGetElement(matrix,2,2) 00196 element(5) = sfGetElement(matrix,2,3) 00197 element(6) = sfGetElement(matrix,3,2) 00198 element(7) = sfGetElement(matrix,3,3) 00199 element(8) = sfGetElement(matrix,3,4) 00200 element(9) = sfGetElement(matrix,4,3) 00201 element(10) = sfGetElement(matrix,4,4) 00202 call sfClear(matrix) 00203 call sfAdd1Real(element(1), 2d0) 00204 call sfAdd1Real(element(2), -1d0) 00205 call sfAdd1Real(element(3), -1d0) 00206 call sfAdd1Real(element(4), 3d0) 00207 call sfAdd1Real(element(5), -1d0) 00208 call sfAdd1Real(element(6), -1d0) 00209 call sfAdd1Real(element(7), 3d0) 00210 call sfAdd1Real(element(8), -1d0) 00211 call sfAdd1Real(element(9), -1d0) 00212 call sfAdd1Real(element(10), 3d0) 00213 call sfprint(matrix, .false., .false.) 00214 rhs(1) = 34d0 00215 rhs(2) = 0d0 00216 rhs(3) = 0d0 00217 rhs(4) = 0d0 00218 error = sfFactor(matrix) 00219 call sfSolve(matrix, rhs, solution) 00220 write (6, 10) rhs(1), rhs(2), rhs(3), rhs(4) 00221 10 format (f 10.2) 00222 end 00223 00224 * 00225 */ 00226 00227 00228 00229 00230 00231 00232 /* 00233 * MATRIX ALLOCATION 00234 * 00235 * Allocates and initializes the data structures associated with a matrix. 00236 * 00237 * >>> Returned: [INTEGER] 00238 * A pointer to the matrix is returned cast into an integer. This pointer 00239 * is then passed and used by the other matrix routines to refer to a 00240 * particular matrix. If an error occurs, the NULL pointer is returned. 00241 * 00242 * >>> Arguments: 00243 * Size <input> (long *) [INTEGER] 00244 * Size of matrix or estimate of size of matrix if matrix is EXPANDABLE. 00245 * Complex <input> (int *) [INTEGER or INTEGER*2] 00246 * Type of matrix. If Complex is 0 then the matrix is real, otherwise 00247 * the matrix will be complex. Note that if the routines are not set up 00248 * to handle the type of matrix requested, then a spPANIC error will occur. 00249 * Further note that if a matrix will be both real and complex, it must 00250 * be specified here as being complex. 00251 * Error <output> (int *) [INTEGER or INTEGER*2] 00252 * Returns error flag, needed because function spError() will not work 00253 * correctly if spCreate() returns NULL. 00254 * 00255 * >>> Possible errors: 00256 * spNO_MEMORY 00257 * spPANIC 00258 * Error is cleared in this routine. 00259 */ 00260 00261 long 00262 sfCreate( Size, Complex, Error ) 00263 00264 int *Size, *Complex, *Error; 00265 { 00266 /* Begin `sfCreate'. */ 00267 return (long)spCreate(*Size, *Complex, Error ); 00268 } 00269 00270 00271 00272 00273 00274 00275 /* 00276 * MATRIX DEALLOCATION 00277 * 00278 * Deallocates pointers and elements of Matrix. 00279 * 00280 * >>> Arguments: 00281 * Matrix <input> (long *) [INTEGER] 00282 * Pointer to the matrix frame which is to be removed from memory. 00283 */ 00284 00285 void 00286 sfDestroy( Matrix ) 00287 00288 long *Matrix; 00289 { 00290 /* Begin `sfDestroy'. */ 00291 spDestroy((char *)*Matrix); 00292 return; 00293 } 00294 00295 00296 00297 00298 00299 00300 #if STRIP 00301 00302 /* 00303 * STRIP FILL-INS FROM MATRIX 00304 * 00305 * Strips the matrix of all fill-ins. 00306 * 00307 * >>> Arguments: 00308 * Matrix <input> (long *) [INTEGER] 00309 * Pointer to the matrix to be stripped. 00310 */ 00311 00312 void 00313 sfStripFills( Matrix ) 00314 00315 long *Matrix; 00316 { 00317 /* Begin `sfStripFills'. */ 00318 spStripFills((char *)*Matrix); 00319 return; 00320 } 00321 #endif 00322 00323 00324 00325 00326 00327 00328 00329 /* 00330 * CLEAR MATRIX 00331 * 00332 * Sets every element of the matrix to zero and clears the error flag. 00333 * 00334 * >>> Arguments: 00335 * Matrix <input> (long *) [INTEGER] 00336 * Pointer to matrix that is to be cleared. 00337 */ 00338 00339 void 00340 sfClear( Matrix ) 00341 00342 long *Matrix; 00343 { 00344 /* Begin `sfClear'. */ 00345 spClear((char *)*Matrix); 00346 return; 00347 } 00348 00349 00350 00351 00352 00353 00354 /* 00355 * SINGLE ELEMENT ADDITION TO MATRIX BY INDEX 00356 * 00357 * Finds element [Row,Col] and returns a pointer to it. If element is 00358 * not found then it is created and spliced into matrix. This routine 00359 * is only to be used after spCreate() and before spMNA_Preorder(), 00360 * spFactor() or spOrderAndFactor(). Returns a pointer to the 00361 * Real portion of a MatrixElement. This pointer is later used by 00362 * sfAddxxxxx() to directly access element. 00363 * 00364 * >>> Returns: [INTEGER] 00365 * Returns a pointer to the element. This pointer is then used to directly 00366 * access the element during successive builds. 00367 * 00368 * >>> Arguments: 00369 * Matrix <input> (long *) [INTEGER] 00370 * Pointer to the matrix that the element is to be added to. 00371 * Row <input> (int *) [INTEGER or INTEGER*2] 00372 * Row index for element. Must be in the range of [0..Size] unless 00373 * the options EXPANDABLE or TRANSLATE are used. Elements placed in 00374 * row zero are discarded. In no case may Row be less than zero. 00375 * Col <input> (int *) [INTEGER or INTEGER*2] 00376 * Column index for element. Must be in the range of [0..Size] unless 00377 * the options EXPANDABLE or TRANSLATE are used. Elements placed in 00378 * column zero are discarded. In no case may Col be less than zero. 00379 * 00380 * >>> Possible errors: 00381 * spNO_MEMORY 00382 * Error is not cleared in this routine. 00383 */ 00384 00385 long 00386 sfGetElement( Matrix, Row, Col ) 00387 00388 long *Matrix; 00389 int *Row, *Col; 00390 { 00391 /* Begin `sfGetElement'. */ 00392 return (long)spGetElement((char *)*Matrix, *Row, *Col); 00393 } 00394 00395 00396 00397 00398 00399 00400 00401 #if QUAD_ELEMENT 00402 /* 00403 * ADDITION OF ADMITTANCE TO MATRIX BY INDEX 00404 * 00405 * Performs same function as sfGetElement except rather than one 00406 * element, all four Matrix elements for a floating component are 00407 * added. This routine also works if component is grounded. Positive 00408 * elements are placed at [Node1,Node2] and [Node2,Node1]. This 00409 * routine is only to be used after sfCreate() and before 00410 * sfMNA_Preorder(), sfFactor() or sfOrderAndFactor(). 00411 * 00412 * >>> Returns: [INTEGER or INTEGER*2] 00413 * Error code. 00414 * 00415 * >>> Arguments: 00416 * Matrix <input> (long *) [INTEGER] 00417 * Pointer to the matrix that component is to be entered in. 00418 * Node1 <input> (int *) [INTEGER or INTEGER*2] 00419 * Row and column indices for elements. Must be in the range of [0..Size] 00420 * unless the options EXPANDABLE or TRANSLATE are used. Node zero is the 00421 * ground node. In no case may Node1 be less than zero. 00422 * Node2 <input> (int *) [INTEGER or INTEGER*2] 00423 * Row and column indices for elements. Must be in the range of [0..Size] 00424 * unless the options EXPANDABLE or TRANSLATE are used. Node zero is the 00425 * ground node. In no case may Node2 be less than zero. 00426 * Template <output> (long[4]) [INTEGER (4)] 00427 * Collection of pointers to four elements that are later used to directly 00428 * address elements. User must supply the template, this routine will 00429 * fill it. 00430 * 00431 * Possible errors: 00432 * spNO_MEMORY 00433 * Error is not cleared in this routine. 00434 */ 00435 00436 int 00437 sfGetAdmittance( Matrix, Node1, Node2, Template ) 00438 00439 long *Matrix, Template[4]; 00440 int *Node1, *Node2; 00441 { 00442 /* Begin `spGetAdmittance'. */ 00443 return 00444 ( spGetAdmittance((char *)*Matrix, *Node1, *Node2, 00445 (struct spTemplate *)Template ) 00446 ); 00447 } 00448 #endif /* QUAD_ELEMENT */ 00449 00450 00451 00452 00453 00454 00455 00456 00457 00458 #if QUAD_ELEMENT 00459 /* 00460 * ADDITION OF FOUR ELEMENTS TO MATRIX BY INDEX 00461 * 00462 * Similar to sfGetAdmittance, except that sfGetAdmittance only 00463 * handles 2-terminal components, whereas sfGetQuad handles simple 00464 * 4-terminals as well. These 4-terminals are simply generalized 00465 * 2-terminals with the option of having the sense terminals different 00466 * from the source and sink terminals. sfGetQuad adds four 00467 * elements to the matrix. Positive elements occur at Row1,Col1 00468 * Row2,Col2 while negative elements occur at Row1,Col2 and Row2,Col1. 00469 * The routine works fine if any of the rows and columns are zero. 00470 * This routine is only to be used after sfCreate() and before 00471 * sfMNA_Preorder(), sfFactor() or sfOrderAndFactor() 00472 * unless TRANSLATE is set true. 00473 * 00474 * >>> Returns: [INTEGER or INTEGER*2] 00475 * Error code. 00476 * 00477 * >>> Arguments: 00478 * Matrix <input> (long *) [INTEGER] 00479 * Pointer to the matrix that component is to be entered in. 00480 * Row1 <input> (int *) [INTEGER or INTEGER*2] 00481 * First row index for elements. Must be in the range of [0..Size] 00482 * unless the options EXPANDABLE or TRANSLATE are used. Zero is the 00483 * ground row. In no case may Row1 be less than zero. 00484 * Row2 <input> (int *) [INTEGER or INTEGER*2] 00485 * Second row index for elements. Must be in the range of [0..Size] 00486 * unless the options EXPANDABLE or TRANSLATE are used. Zero is the 00487 * ground row. In no case may Row2 be less than zero. 00488 * Col1 <input> (int *) [INTEGER or INTEGER*2] 00489 * First column index for elements. Must be in the range of [0..Size] 00490 * unless the options EXPANDABLE or TRANSLATE are used. Zero is the 00491 * ground column. In no case may Col1 be less than zero. 00492 * Col2 <input> (int *) [INTEGER or INTEGER*2] 00493 * Second column index for elements. Must be in the range of [0..Size] 00494 * unless the options EXPANDABLE or TRANSLATE are used. Zero is the 00495 * ground column. In no case may Col2 be less than zero. 00496 * Template <output> (long[4]) [INTEGER (4)] 00497 * Collection of pointers to four elements that are later used to directly 00498 * address elements. User must supply the template, this routine will 00499 * fill it. 00500 * 00501 * Possible errors: 00502 * spNO_MEMORY 00503 * Error is not cleared in this routine. 00504 */ 00505 00506 int 00507 sfGetQuad( Matrix, Row1, Row2, Col1, Col2, Template ) 00508 00509 long *Matrix, Template[4]; 00510 int *Row1, *Row2, *Col1, *Col2; 00511 { 00512 /* Begin `spGetQuad'. */ 00513 return 00514 ( spGetQuad( (char *)*Matrix, *Row1, *Row2, *Col1, *Col2, 00515 (struct spTemplate *)Template ) 00516 ); 00517 } 00518 #endif /* QUAD_ELEMENT */ 00519 00520 00521 00522 00523 00524 00525 00526 00527 00528 #if QUAD_ELEMENT 00529 /* 00530 * ADDITION OF FOUR STRUCTURAL ONES TO MATRIX BY INDEX 00531 * 00532 * Performs similar function to sfGetQuad() except this routine is 00533 * meant for components that do not have an admittance representation. 00534 * 00535 * The following stamp is used: 00536 * Pos Neg Eqn 00537 * Pos [ . . 1 ] 00538 * Neg [ . . -1 ] 00539 * Eqn [ 1 -1 . ] 00540 * 00541 * >>> Returns: [INTEGER or INTEGER*2] 00542 * Error code. 00543 * 00544 * >>> Arguments: 00545 * Matrix <input> (long *) [INTEGER] 00546 * Pointer to the matrix that component is to be entered in. 00547 * Pos <input> (int *) [INTEGER or INTEGER*2] 00548 * See stamp above. Must be in the range of [0..Size] 00549 * unless the options EXPANDABLE or TRANSLATE are used. Zero is the 00550 * ground row. In no case may Pos be less than zero. 00551 * Neg <input> (int *) [INTEGER or INTEGER*2] 00552 * See stamp above. Must be in the range of [0..Size] 00553 * unless the options EXPANDABLE or TRANSLATE are used. Zero is the 00554 * ground row. In no case may Neg be less than zero. 00555 * Eqn <input> (int *) [INTEGER or INTEGER*2] 00556 * See stamp above. Must be in the range of [0..Size] 00557 * unless the options EXPANDABLE or TRANSLATE are used. Zero is the 00558 * ground row. In no case may Eqn be less than zero. 00559 * Template <output> (long[4]) [INTEGER (4)] 00560 * Collection of pointers to four elements that are later used to directly 00561 * address elements. User must supply the template, this routine will 00562 * fill it. 00563 * 00564 * Possible errors: 00565 * spNO_MEMORY 00566 * Error is not cleared in this routine. 00567 */ 00568 00569 int 00570 sfGetOnes(Matrix, Pos, Neg, Eqn, Template) 00571 00572 long *Matrix, Template[4]; 00573 int *Pos, *Neg, *Eqn; 00574 { 00575 /* Begin `sfGetOnes'. */ 00576 return 00577 ( spGetOnes( (char *)*Matrix, *Pos, *Neg, *Eqn, 00578 (struct spTemplate *)Template ) 00579 ); 00580 } 00581 #endif /* QUAD_ELEMENT */ 00582 00583 00584 00585 00586 00587 00588 00589 /* 00590 * ADD ELEMENT(S) DIRECTLY TO MATRIX 00591 * 00592 * Adds a value to an element or a set of four element in a matrix. 00593 * These elements are referenced by pointer, and so must already have 00594 * been created by spGetElement(), spGetAdmittance(), spGetQuad(), or 00595 * spGetOnes(). 00596 * 00597 * >>> Arguments: 00598 * Element <input> (long *) [INTEGER] 00599 * Pointer to the element that is to be added to. 00600 * Template <input> (long[4]) [INTEGER (4)] 00601 * Pointer to the element that is to be added to. 00602 * Real <input> (double *) [REAL or DOUBLE PRECISION] 00603 * Real portion of the number to be added to the element. 00604 * Imag <input> (double *) [REAL or DOUBLE PRECISION] 00605 * Imaginary portion of the number to be added to the element. 00606 */ 00607 00608 void 00609 sfAdd1Real( Element, Real ) 00610 00611 long *Element; 00612 RealNumber *Real; 00613 { 00614 /* Begin `sfAdd1Real'. */ 00615 *((RealNumber *)*Element) += *Real; 00616 } 00617 00618 00619 #if spCOMPLEX 00620 00621 void 00622 sfAdd1Imag( Element, Imag ) 00623 00624 long *Element; 00625 RealNumber *Imag; 00626 { 00627 /* Begin `sfAdd1Imag'. */ 00628 *(((RealNumber *)*Element)+1) += *Imag; 00629 } 00630 00631 00632 void 00633 sfAdd1Complex( Element, Real, Imag ) 00634 00635 long *Element; 00636 RealNumber *Real, *Imag; 00637 { 00638 /* Begin `sfAdd1Complex'. */ 00639 *((RealNumber *)*Element) += *Real; 00640 *(((RealNumber *)*Element)+1) += *Imag; 00641 } 00642 #endif /* spCOMPLEX */ 00643 00644 00645 #if QUAD_ELEMENT 00646 00647 void 00648 sfAdd4Real( Template, Real ) 00649 00650 long Template[4]; 00651 RealNumber *Real; 00652 { 00653 /* Begin `sfAdd4Real'. */ 00654 *((RealNumber *)Template[0]) += *Real; 00655 *((RealNumber *)Template[1]) += *Real; 00656 *((RealNumber *)Template[2]) -= *Real; 00657 *((RealNumber *)Template[3]) -= *Real; 00658 } 00659 00660 00661 #if spCOMPLEX 00662 00663 void 00664 sfAdd4Imag( Template, Imag ) 00665 00666 long Template[4]; 00667 RealNumber *Imag; 00668 { 00669 /* Begin `sfAdd4Imag'. */ 00670 *(((RealNumber *)Template[0])+1) += *Imag; 00671 *(((RealNumber *)Template[1])+1) += *Imag; 00672 *(((RealNumber *)Template[2])+1) -= *Imag; 00673 *(((RealNumber *)Template[3])+1) -= *Imag; 00674 } 00675 00676 00677 void 00678 sfAdd4Complex( Template, Real, Imag ) 00679 00680 long Template[4]; 00681 RealNumber *Real, *Imag; 00682 { 00683 /* Begin `sfAdd4Complex'. */ 00684 *((RealNumber *)Template[0]) += *Real; 00685 *((RealNumber *)Template[1]) += *Real; 00686 *((RealNumber *)Template[2]) -= *Real; 00687 *((RealNumber *)Template[3]) -= *Real; 00688 *(((RealNumber *)Template[0])+1) += *Imag; 00689 *(((RealNumber *)Template[1])+1) += *Imag; 00690 *(((RealNumber *)Template[2])+1) -= *Imag; 00691 *(((RealNumber *)Template[3])+1) -= *Imag; 00692 } 00693 #endif /* spCOMPLEX */ 00694 #endif /* QUAD_ELEMENT */ 00695 00696 00697 00698 00699 00700 00701 /* 00702 * ORDER && FACTOR MATRIX 00703 * 00704 * This routine chooses a pivot order for the matrix and factors it 00705 * into LU form. It handles both the initial factorization and subsequent 00706 * factorizations when a reordering is desired. This is handled in a manner 00707 * that is transparent to the user. The routine uses a variation of 00708 * Gauss's method where the pivots are associated with L and the 00709 * diagonal terms of U are one. 00710 * 00711 * >>> Returned: [INTEGER of INTEGER*2] 00712 * The error code is returned. Possible errors are listed below. 00713 * 00714 * >>> Arguments: 00715 * Matrix <input> (long *) [INTEGER] 00716 * Pointer to matrix. 00717 * RHS <input> (RealVector) [REAL (1) or DOUBLE PRECISION (1)] 00718 * Representative right-hand side vector that is used to determine 00719 * pivoting order when the right hand side vector is sparse. If 00720 * RHS is a NULL pointer then the RHS vector is assumed to 00721 * be full and it is not used when determining the pivoting 00722 * order. 00723 * RelThreshold <input> (RealNumber *) [REAL or DOUBLE PRECISION] 00724 * This number determines what the pivot relative threshold will 00725 * be. It should be between zero and one. If it is one then the 00726 * pivoting method becomes complete pivoting, which is very slow 00727 * and tends to fill up the matrix. If it is set close to zero 00728 * the pivoting method becomes strict Markowitz with no 00729 * threshold. The pivot threshold is used to eliminate pivot 00730 * candidates that would cause excessive element growth if they 00731 * were used. Element growth is the cause of roundoff error. 00732 * Element growth occurs even in well-conditioned matrices. 00733 * Setting the RelThreshold large will reduce element growth and 00734 * roundoff error, but setting it too large will cause execution 00735 * time to be excessive and will result in a large number of 00736 * fill-ins. If this occurs, accuracy can actually be degraded 00737 * because of the large number of operations required on the 00738 * matrix due to the large number of fill-ins. A good value seems 00739 * to be 0.001. The default is chosen by giving a value larger 00740 * than one or less than or equal to zero. This value should be 00741 * increased and the matrix resolved if growth is found to be 00742 * excessive. Changing the pivot threshold does not improve 00743 * performance on matrices where growth is low, as is often the 00744 * case with ill-conditioned matrices. Once a valid threshold is 00745 * given, it becomes the new default. The default value of 00746 * RelThreshold was choosen for use with nearly diagonally 00747 * dominant matrices such as node- and modified-node admittance 00748 * matrices. For these matrices it is usually best to use 00749 * diagonal pivoting. For matrices without a strong diagonal, it 00750 * is usually best to use a larger threshold, such as 0.01 or 00751 * 0.1. 00752 * AbsThreshold <input> (RealNumber *) [REAL or DOUBLE PRECISION] 00753 * The absolute magnitude an element must have to be considered 00754 * as a pivot candidate, except as a last resort. This number 00755 * should be set significantly smaller than the smallest diagonal 00756 * element that is is expected to be placed in the matrix. If 00757 * there is no reasonable prediction for the lower bound on these 00758 * elements, then AbsThreshold should be set to zero. 00759 * AbsThreshold is used to reduce the possibility of choosing as a 00760 * pivot an element that has suffered heavy cancellation and as a 00761 * result mainly consists of roundoff error. Once a valid 00762 * threshold is given, it becomes the new default. 00763 * DiagPivoting <input> (long *) [LOGICAL] 00764 * A flag indicating that pivot selection should be confined to the 00765 * diagonal if possible. If DiagPivoting is nonzero and if 00766 * DIAGONAL_PIVOTING is enabled pivots will be chosen only from 00767 * the diagonal unless there are no diagonal elements that satisfy 00768 * the threshold criteria. Otherwise, the entire reduced 00769 * submatrix is searched when looking for a pivot. The diagonal 00770 * pivoting in Sparse is efficient and well refined, while the 00771 * off-diagonal pivoting is not. For symmetric and near symmetric 00772 * matrices, it is best to use diagonal pivoting because it 00773 * results in the best performance when reordering the matrix and 00774 * when factoring the matrix without ordering. If there is a 00775 * considerable amount of nonsymmetry in the matrix, then 00776 * off-diagonal pivoting may result in a better equation ordering 00777 * simply because there are more pivot candidates to choose from. 00778 * A better ordering results in faster subsequent factorizations. 00779 * However, the initial pivot selection process takes considerably 00780 * longer for off-diagonal pivoting. 00781 * 00782 * >>> Possible errors: 00783 * spNO_MEMORY 00784 * spSINGULAR 00785 * spSMALL_PIVOT 00786 * Error is cleared in this function. 00787 */ 00788 00789 int 00790 sfOrderAndFactor( Matrix, RHS, RelThreshold, AbsThreshold, DiagPivoting ) 00791 00792 long *Matrix, *DiagPivoting; 00793 RealNumber RHS[], *RelThreshold, *AbsThreshold; 00794 { 00795 /* Begin `sfOrderAndFactor'. */ 00796 return spOrderAndFactor( (char *)*Matrix, RHS, *RelThreshold, 00797 *AbsThreshold, (BOOLEAN)*DiagPivoting ); 00798 } 00799 00800 00801 00802 00803 00804 00805 00806 /* 00807 * FACTOR MATRIX 00808 * 00809 * This routine is the companion routine to spOrderAndFactor(). 00810 * Unlike sfOrderAndFactor(), sfFactor() cannot change the ordering. 00811 * It is also faster than sfOrderAndFactor(). The standard way of 00812 * using these two routines is to first use sfOrderAndFactor() for the 00813 * initial factorization. For subsequent factorizations, sfFactor() 00814 * is used if there is some assurance that little growth will occur 00815 * (say for example, that the matrix is diagonally dominant). If 00816 * sfFactor() is called for the initial factorization of the matrix, 00817 * then sfOrderAndFactor() is automatically called with the default 00818 * threshold. This routine uses "row at a time" LU factorization. 00819 * Pivots are associated with the lower triangular matrix and the 00820 * diagonals of the upper triangular matrix are ones. 00821 * 00822 * >>> Returned: [INTEGER or INTEGER*2] 00823 * The error code is returned. Possible errors are listed below. 00824 * 00825 * >>> Arguments: 00826 * Matrix <input> (long *) [INTEGER] 00827 * Pointer to matrix. 00828 * 00829 * >>> Possible errors: 00830 * spNO_MEMORY 00831 * spSINGULAR 00832 * spZERO_DIAG 00833 * spSMALL_PIVOT 00834 * Error is cleared in this function. 00835 */ 00836 00837 int 00838 sfFactor( Matrix ) 00839 00840 long *Matrix; 00841 { 00842 /* Begin `sfFactor'. */ 00843 return spFactor((char *)*Matrix); 00844 } 00845 00846 00847 00848 00849 00850 00851 /* 00852 * PARTITION MATRIX 00853 * 00854 * This routine determines the cost to factor each row using both 00855 * direct and indirect addressing and decides, on a row-by-row basis, 00856 * which addressing mode is fastest. This information is used in 00857 * sfFactor() to speed the factorization. 00858 * 00859 * When factoring a previously ordered matrix using sfFactor(), Sparse 00860 * operates on a row-at-a-time basis. For speed, on each step, the 00861 * row being updated is copied into a full vector and the operations 00862 * are performed on that vector. This can be done one of two ways, 00863 * either using direct addressing or indirect addressing. Direct 00864 * addressing is fastest when the matrix is relatively dense and 00865 * indirect addressing is best when the matrix is quite sparse. The 00866 * user selects the type of partition used with Mode. If Mode is set 00867 * to spDIRECT_PARTITION, then the all rows are placed in the direct 00868 * addressing partition. Similarly, if Mode is set to 00869 * spINDIRECT_PARTITION, then the all rows are placed in the indirect 00870 * addressing partition. By setting Mode to spAUTO_PARTITION, the 00871 * user allows Sparse to select the partition for each row 00872 * individually. sfFactor() generally runs faster if Sparse is 00873 * allowed to choose its own partitioning, however choosing a 00874 * partition is expensive. The time required to choose a partition is 00875 * of the same order of the cost to factor the matrix. If you plan to 00876 * factor a large number of matrices with the same structure, it is 00877 * best to let Sparse choose the partition. Otherwise, you should 00878 * choose the partition based on the predicted density of the matrix. 00879 * 00880 * >>> Arguments: 00881 * Matrix <input> (long *) [INTEGER] 00882 * Pointer to matrix. 00883 * Mode <input> (int *) [INTEGER or INTEGER*2] 00884 * Mode must be one of three special codes: spDIRECT_PARTITION, 00885 * spINDIRECT_PARTITION, or spAUTO_PARTITION. 00886 */ 00887 00888 void 00889 sfPartition( Matrix, Mode ) 00890 00891 long *Matrix; 00892 int *Mode; 00893 { 00894 /* Begin `sfPartition'. */ 00895 spPartition((char *)*Matrix, *Mode); 00896 } 00897 00898 00899 00900 00901 00902 00903 00904 /* 00905 * SOLVE MATRIX EQUATION 00906 * 00907 * Performs forward elimination and back substitution to find the 00908 * unknown vector from the RHS vector and factored matrix. This 00909 * routine assumes that the pivots are associated with the lower 00910 * triangular (L) matrix and that the diagonal of the upper triangular 00911 * (U) matrix consists of ones. This routine arranges the computation 00912 * in different way than is traditionally used in order to exploit the 00913 * sparsity of the right-hand side. See the reference in spRevision. 00914 * 00915 * >>> Arguments: 00916 * Matrix <input> (long *) [INTEGER] 00917 * Pointer to matrix. 00918 * RHS <input> (RealVector) [REAL (1) or DOUBLE PRECISION (1)] 00919 * RHS is the input data array, the right hand side. This data is 00920 * undisturbed and may be reused for other solves. 00921 * Solution <output> (RealVector) [REAL (1) or DOUBLE PRECISION (1)] 00922 * Solution is the output data array. This routine is constructed such that 00923 * RHS and Solution can be the same array. 00924 * iRHS <input> (RealVector) [REAL (1) or DOUBLE PRECISION (1)] 00925 * iRHS is the imaginary portion of the input data array, the right 00926 * hand side. This data is undisturbed and may be reused for other solves. 00927 * This argument is only necessary if matrix is complex and if 00928 * spSEPARATED_COMPLEX_VECTOR is set true. 00929 * iSolution <output> (RealVector) [REAL (1) or DOUBLE PRECISION (1)] 00930 * iSolution is the imaginary portion of the output data array. This 00931 * routine is constructed such that iRHS and iSolution can be 00932 * the same array. This argument is only necessary if matrix is complex 00933 * and if spSEPARATED_COMPLEX_VECTOR is set true. 00934 * 00935 * >>> Obscure Macros 00936 * IMAG_VECTORS 00937 * Replaces itself with `, iRHS, iSolution' if the options spCOMPLEX and 00938 * spSEPARATED_COMPLEX_VECTORS are set, otherwise it disappears 00939 * without a trace. 00940 */ 00941 00942 /*VARARGS3*/ 00943 00944 void 00945 sfSolve( Matrix, RHS, Solution IMAG_VECTORS ) 00946 00947 long *Matrix; 00948 RealVector RHS, Solution IMAG_VECTORS; 00949 { 00950 /* Begin `sfSolve'. */ 00951 spSolve( (char *)*Matrix, RHS, Solution IMAG_VECTORS ); 00952 } 00953 00954 00955 00956 00957 00958 00959 #if TRANSPOSE 00960 /* 00961 * SOLVE TRANSPOSED MATRIX EQUATION 00962 * 00963 * Performs forward elimination and back substitution to find the 00964 * unknown vector from the RHS vector and transposed factored 00965 * matrix. This routine is useful when performing sensitivity analysis 00966 * on a circuit using the adjoint method. This routine assumes that 00967 * the pivots are associated with the untransposed lower triangular 00968 * (L) matrix and that the diagonal of the untransposed upper 00969 * triangular (U) matrix consists of ones. 00970 * 00971 * >>> Arguments: 00972 * Matrix <input> (long *) [INTEGER] 00973 * Pointer to matrix. 00974 * RHS <input> (RealVector) [REAL (1) or DOUBLE PRECISION (1)] 00975 * RHS is the input data array, the right hand side. This data is 00976 * undisturbed and may be reused for other solves. 00977 * Solution <output> (RealVector) [REAL (1) or DOUBLE PRECISION (1)] 00978 * Solution is the output data array. This routine is constructed such that 00979 * RHS and Solution can be the same array. 00980 * iRHS <input> (RealVector) [REAL (1) or DOUBLE PRECISION (1)] 00981 * iRHS is the imaginary portion of the input data array, the right 00982 * hand side. This data is undisturbed and may be reused for other solves. 00983 * If spSEPARATED_COMPLEX_VECTOR is set false, or if matrix is real, there 00984 * is no need to supply this array. 00985 * iSolution <output> (RealVector) [REAL (1) or DOUBLE PRECISION (1)] 00986 * iSolution is the imaginary portion of the output data array. This 00987 * routine is constructed such that iRHS and iSolution can be 00988 * the same array. If spSEPARATED_COMPLEX_VECTOR is set false, or if 00989 * matrix is real, there is no need to supply this array. 00990 * 00991 * >>> Obscure Macros 00992 * IMAG_VECTORS 00993 * Replaces itself with `, iRHS, iSolution' if the options spCOMPLEX and 00994 * spSEPARATED_COMPLEX_VECTORS are set, otherwise it disappears 00995 * without a trace. 00996 */ 00997 00998 /*VARARGS3*/ 00999 01000 void 01001 sfSolveTransposed( Matrix, RHS, Solution IMAG_VECTORS ) 01002 01003 long *Matrix; 01004 RealVector RHS, Solution IMAG_VECTORS; 01005 { 01006 /* Begin `sfSolveTransposed'. */ 01007 spSolveTransposed( (char *)*Matrix, RHS, Solution IMAG_VECTORS ); 01008 } 01009 #endif /* TRANSPOSE */ 01010 01011 01012 01013 01014 01015 #if DOCUMENTATION 01016 /* 01017 * PRINT MATRIX 01018 * 01019 * Formats and send the matrix to standard output. Some elementary 01020 * statistics are also output. The matrix is output in a format that is 01021 * readable by people. 01022 * 01023 * >>> Arguments: 01024 * Matrix <input> (long *) [INTEGER] 01025 * Pointer to matrix. 01026 * PrintReordered <input> (long *) [LOGICAL] 01027 * Indicates whether the matrix should be printed out in its original 01028 * form, as input by the user, or whether it should be printed in its 01029 * reordered form, as used by the matrix routines. A zero indicates that 01030 * the matrix should be printed as inputed, a one indicates that it 01031 * should be printed reordered. 01032 * Data <input> (long *) [LOGICAL] 01033 * Boolean flag that when false indicates that output should be 01034 * compressed such that only the existence of an element should be 01035 * indicated rather than giving the actual value. Thus 10 times as many 01036 * can be printed on a row. A zero signifies that the matrix should 01037 * be printed compressed. A one indicates that the matrix should be 01038 * printed in all its glory. 01039 * Header <input> (long *) [LOGICAL] 01040 * Flag indicating that extra information such as the row and column 01041 * numbers should be printed. 01042 */ 01043 01044 void 01045 sfPrint( Matrix, Data, PrintReordered, Header ) 01046 01047 long *Matrix, *PrintReordered, *Data, *Header; 01048 { 01049 /* Begin `sfPrint'. */ 01050 spPrint( (char *)*Matrix, (int)*PrintReordered, (int)*Data, (int)*Header ); 01051 } 01052 #endif /* DOCUMENTATION */ 01053 01054 01055 01056 01057 01058 01059 #if DOCUMENTATION 01060 /* 01061 * OUTPUT MATRIX TO FILE 01062 * 01063 * Writes matrix to file in format suitable to be read back in by the 01064 * matrix test program. Data is sent to a file with a fixed name because 01065 * it is impossible to pass strings from FORTRAN to C in a manner that is 01066 * portable. 01067 * 01068 * >>> Returns: 01069 * One is returned if routine was successful, otherwise zero is returned. 01070 * The calling function can query errno (the system global error variable) 01071 * as to the reason why this routine failed. 01072 * 01073 * >>> Arguments: [LOGICAL] 01074 * Matrix <input> (long *) [INTEGER] 01075 * Pointer to matrix. 01076 * Reordered <input> (long *) [LOGICAL] 01077 * Specifies whether matrix should be output in reordered form, 01078 * or in original order. 01079 * Data <input> (long *) [LOGICAL] 01080 * Indicates that the element values should be output along with 01081 * the indices for each element. This parameter must be true if 01082 * matrix is to be read by the sparse test program. 01083 * Header <input> (long *) [LOGICAL] 01084 * Indicates that header is desired. This parameter must be true if 01085 * matrix is to be read by the sparse test program. 01086 */ 01087 #define MATRIX_FILE_NAME "spMatrix" 01088 #define STATS_FILE_NAME "spStats" 01089 01090 long 01091 sfFileMatrix( Matrix, Reordered, Data, Header ) 01092 01093 long *Matrix, *Reordered, *Data, *Header; 01094 { 01095 /* Begin `sfFileMatrix'. */ 01096 return spFileMatrix( (char *)*Matrix, MATRIX_FILE_NAME, "", 01097 (int)*Reordered, (int)*Data, (int)*Header ); 01098 } 01099 #endif /* DOCUMENTATION */ 01100 01101 01102 01103 01104 01105 01106 #if DOCUMENTATION 01107 /* 01108 * OUTPUT SOURCE VECTOR TO FILE 01109 * 01110 * Writes vector to file in format suitable to be read back in by the 01111 * matrix test program. This routine should be executed after the function 01112 * sfFileMatrix. 01113 * 01114 * >>> Returns: 01115 * One is returned if routine was successful, otherwise zero is returned. 01116 * The calling function can query errno (the system global error variable) 01117 * as to the reason why this routine failed. 01118 * 01119 * >>> Arguments: 01120 * Matrix <input> (long *) 01121 * Pointer to matrix. 01122 * RHS <input> (RealNumber []) [REAL (1) or DOUBLE PRECISION (1)] 01123 * Right-hand side vector. This is only the real portion if 01124 * spSEPARATED_COMPLEX_VECTORS is true. 01125 * iRHS <input> (RealNumber []) [REAL (1) or DOUBLE PRECISION (1)] 01126 * Right-hand side vector, imaginary portion. Not necessary if matrix 01127 * is real or if spSEPARATED_COMPLEX_VECTORS is set false. 01128 */ 01129 01130 int 01131 sfFileVector( Matrix, RHS IMAG_RHS ) 01132 01133 long *Matrix; 01134 RealVector RHS IMAG_RHS; 01135 { 01136 /* Begin `sfFileVector'. */ 01137 return spFileVector( (char *)*Matrix, MATRIX_FILE_NAME, RHS IMAG_RHS ); 01138 } 01139 #endif /* DOCUMENTATION */ 01140 01141 01142 01143 01144 01145 01146 01147 #if DOCUMENTATION 01148 /* 01149 * OUTPUT STATISTICS TO FILE 01150 * 01151 * Writes useful information concerning the matrix to a file. Should be 01152 * executed after the matrix is factored. 01153 * 01154 * >>> Returns: [LOGICAL] 01155 * One is returned if routine was successful, otherwise zero is returned. 01156 * The calling function can query errno (the system global error variable) 01157 * as to the reason why this routine failed. 01158 * 01159 * >>> Arguments: 01160 * Matrix <input> (long *) [INTEGER] 01161 * Pointer to matrix. 01162 */ 01163 01164 int 01165 sfFileStats( Matrix ) 01166 01167 long *Matrix; 01168 { 01169 /* Begin `sfFileStats'. */ 01170 return spFileStats( (char *)*Matrix, STATS_FILE_NAME, "" ); 01171 } 01172 #endif /* DOCUMENTATION */ 01173 01174 01175 01176 01177 #if MODIFIED_NODAL 01178 /* 01179 * PREORDER MODIFIED NODE ADMITTANCE MATRIX TO REMOVE ZEROS FROM DIAGONAL 01180 * 01181 * This routine massages modified node admittance matrices to remove 01182 * zeros from the diagonal. It takes advantage of the fact that the 01183 * row and column associated with a zero diagonal usually have 01184 * structural ones placed symmetricly. This routine should be used 01185 * only on modified node admittance matrices and should be executed 01186 * after the matrix has been built but before the factorization 01187 * begins. It should be executed for the initial factorization only 01188 * and should be executed before the rows have been linked. Thus it 01189 * should be run before using spScale(), spMultiply(), 01190 * spDeleteRowAndCol(), or spNorm(). 01191 * 01192 * This routine exploits the fact that the structural one are placed 01193 * in the matrix in symmetric twins. For example, the stamps for 01194 * grounded and a floating voltage sources are 01195 * grounded: floating: 01196 * [ x x 1 ] [ x x 1 ] 01197 * [ x x ] [ x x -1 ] 01198 * [ 1 ] [ 1 -1 ] 01199 * Notice for the grounded source, there is one set of twins, and for 01200 * the grounded, there are two sets. We remove the zero from the diagonal 01201 * by swapping the rows associated with a set of twins. For example: 01202 * grounded: floating 1: floating 2: 01203 * [ 1 ] [ 1 -1 ] [ x x 1 ] 01204 * [ x x ] [ x x -1 ] [ 1 -1 ] 01205 * [ x x 1 ] [ x x 1 ] [ x x -1 ] 01206 * 01207 * It is important to deal with any zero diagonals that only have one 01208 * set of twins before dealing with those that have more than one because 01209 * swapping row destroys the symmetry of any twins in the rows being 01210 * swapped, which may limit future moves. Consider 01211 * [ x x 1 ] 01212 * [ x x -1 1 ] 01213 * [ 1 -1 ] 01214 * [ 1 ] 01215 * There is one set of twins for diagonal 4 and two for diagonal3. 01216 * Dealing with diagonal for first requires swapping rows 2 and 4. 01217 * [ x x 1 ] 01218 * [ 1 ] 01219 * [ 1 -1 ] 01220 * [ x x -1 1 ] 01221 * We can now deal with diagonal 3 by swapping rows 1 and 3. 01222 * [ 1 -1 ] 01223 * [ 1 ] 01224 * [ x x 1 ] 01225 * [ x x -1 1 ] 01226 * And we are done, there are no zeros left on the diagonal. However, if 01227 * we originally dealt with diagonal 3 first, we could swap rows 2 and 3 01228 * [ x x 1 ] 01229 * [ 1 -1 ] 01230 * [ x x -1 1 ] 01231 * [ 1 ] 01232 * Diagonal 4 no longer has a symmetric twin and we cannot continue. 01233 * 01234 * So we always take care of lone twins first. When none remain, we 01235 * choose arbitrarily a set of twins for a diagonal with more than one set 01236 * and swap the rows corresponding to that twin. We then deal with any 01237 * lone twins that were created and repeat the procedure until no 01238 * zero diagonals with symmetric twins remain. 01239 * 01240 * In this particular implementation, columns are swapped rather than rows. 01241 * The algorithm used in this function was developed by Ken Kundert and 01242 * Tom Quarles. 01243 * 01244 * >>> Arguments: 01245 * Matrix <input> (long *) [INTEGER] 01246 * Pointer to the matrix to be preordered. 01247 */ 01248 01249 void 01250 sfMNA_Preorder( Matrix ) 01251 01252 long *Matrix; 01253 { 01254 /* Begin `sfMNA_Preorder'. */ 01255 spMNA_Preorder( (char *)*Matrix ); 01256 } 01257 #endif /* MODIFIED_NODAL */ 01258 01259 01260 01261 01262 01263 01264 #if SCALING 01265 /* 01266 * SCALE MATRIX 01267 * 01268 * This function scales the matrix to enhance the possibility of 01269 * finding a good pivoting order. Note that scaling enhances accuracy 01270 * of the solution only if it affects the pivoting order, so it makes 01271 * no sense to scale the matrix before spFactor(). If scaling is 01272 * desired it should be done before spOrderAndFactor(). There 01273 * are several things to take into account when choosing the scale 01274 * factors. First, the scale factors are directly multiplied against 01275 * the elements in the matrix. To prevent roundoff, each scale factor 01276 * should be equal to an integer power of the number base of the 01277 * machine. Since most machines operate in base two, scale factors 01278 * should be a power of two. Second, the matrix should be scaled such 01279 * that the matrix of element uncertainties is equilibrated. Third, 01280 * this function multiplies the scale factors by the elements, so if 01281 * one row tends to have uncertainties 1000 times smaller than the 01282 * other rows, then its scale factor should be 1024, not 1/1024. 01283 * Fourth, to save time, this function does not scale rows or columns 01284 * if their scale factors are equal to one. Thus, the scale factors 01285 * should be normalized to the most common scale factor. Rows and 01286 * columns should be normalized separately. For example, if the size 01287 * of the matrix is 100 and 10 rows tend to have uncertainties near 01288 * 1e-6 and the remaining 90 have uncertainties near 1e-12, then the 01289 * scale factor for the 10 should be 1/1,048,576 and the scale factors 01290 * for the remaining 90 should be 1. Fifth, since this routine 01291 * directly operates on the matrix, it is necessary to apply the scale 01292 * factors to the RHS and Solution vectors. It may be easier to 01293 * simply use spOrderAndFactor() on a scaled matrix to choose the 01294 * pivoting order, and then throw away the matrix. Subsequent 01295 * factorizations, performed with spFactor(), will not need to have 01296 * the RHS and Solution vectors descaled. Lastly, this function 01297 * should not be executed before the function spMNA_Preorder. 01298 * 01299 * >>> Arguments: 01300 * Matrix <input> (long *) [INTEGER] 01301 * Pointer to the matrix to be scaled. 01302 * SolutionScaleFactors <input> (RealVector) [REAL(1) or DOUBLE PRECISION(1)] 01303 * The array of Solution scale factors. These factors scale the columns. 01304 * All scale factors are real valued. 01305 * RHS_ScaleFactors <input> (RealVector) [REAL(1) or DOUBLE PRECISION(1)] 01306 * The array of RHS scale factors. These factors scale the rows. 01307 * All scale factors are real valued. 01308 */ 01309 01310 void 01311 sfScale( Matrix, RHS_ScaleFactors, SolutionScaleFactors ) 01312 01313 long *Matrix; 01314 RealVector RHS_ScaleFactors, SolutionScaleFactors; 01315 { 01316 /* Begin `sfScale'. */ 01317 spScale( (char *)*Matrix, RHS_ScaleFactors, SolutionScaleFactors ); 01318 } 01319 #endif /* SCALING */ 01320 01321 01322 01323 01324 01325 01326 #if MULTIPLICATION 01327 /* 01328 * MATRIX MULTIPLICATION 01329 * 01330 * Multiplies matrix by solution vector to find source vector. 01331 * Assumes matrix has not been factored. This routine can be used 01332 * as a test to see if solutions are correct. It should not be used 01333 * before PreorderFoModifiedNodal(). 01334 * 01335 * >>> Arguments: 01336 * Matrix <input> (long *) [INTEGER] 01337 * Pointer to the matrix. 01338 * RHS <output> (RealVector) [REAL(1) or DOUBLE PRECISION(1)] 01339 * RHS is the right hand side. This is what is being solved for. 01340 * Solution <input> (RealVector) [REAL(1) or DOUBLE PRECISION(1)] 01341 * Solution is the vector being multiplied by the matrix. 01342 * iRHS <output> (RealVector) [REAL(1) or DOUBLE PRECISION(1)] 01343 * iRHS is the imaginary portion of the right hand side. This is 01344 * what is being solved for. This is only necessary if the matrix is 01345 * complex and spSEPARATED_COMPLEX_VECTORS is true. 01346 * iSolution <input> (RealVector) [REAL(1) or DOUBLE PRECISION(1)] 01347 * iSolution is the imaginary portion of the vector being multiplied 01348 * by the matrix. This is only necessary if the matrix is 01349 * complex and spSEPARATED_COMPLEX_VECTORS is true. 01350 * 01351 * >>> Obscure Macros 01352 * IMAG_VECTORS 01353 * Replaces itself with `, iRHS, iSolution' if the options spCOMPLEX and 01354 * spSEPARATED_COMPLEX_VECTORS are set, otherwise it disappears 01355 * without a trace. 01356 */ 01357 01358 void 01359 sfMultiply( Matrix, RHS, Solution IMAG_VECTORS ) 01360 01361 long *Matrix; 01362 RealVector Solution, RHS IMAG_VECTORS; 01363 { 01364 /* Begin `sfMultiply'. */ 01365 spMultiply( (char *)*Matrix, RHS, Solution IMAG_VECTORS ); 01366 } 01367 #endif /* MULTIPLICATION */ 01368 01369 01370 01371 01372 01373 01374 #if MULTIPLICATION && TRANSPOSE 01375 /* 01376 * TRANSPOSED MATRIX MULTIPLICATION 01377 * 01378 * Multiplies transposed matrix by solution vector to find source vector. 01379 * Assumes matrix has not been factored. This routine can be used 01380 * as a test to see if solutions are correct. It should not be used 01381 * before PreorderFoModifiedNodal(). 01382 * 01383 * >>> Arguments: 01384 * Matrix <input> (long *) [INTEGER] 01385 * Pointer to the matrix. 01386 * RHS <output> (RealVector) [REAL(1) or DOUBLE PRECISION(1)] 01387 * RHS is the right hand side. This is what is being solved for. 01388 * Solution <input> (RealVector) [REAL(1) or DOUBLE PRECISION(1)] 01389 * Solution is the vector being multiplied by the matrix. 01390 * iRHS <output> (RealVector) [REAL(1) or DOUBLE PRECISION(1)] 01391 * iRHS is the imaginary portion of the right hand side. This is 01392 * what is being solved for. This is only necessary if the matrix is 01393 * complex and spSEPARATED_COMPLEX_VECTORS is true. 01394 * iSolution <input> (RealVector) [REAL(1) or DOUBLE PRECISION(1)] 01395 * iSolution is the imaginary portion of the vector being multiplied 01396 * by the matrix. This is only necessary if the matrix is 01397 * complex and spSEPARATED_COMPLEX_VECTORS is true. 01398 * 01399 * >>> Obscure Macros 01400 * IMAG_VECTORS 01401 * Replaces itself with `, iRHS, iSolution' if the options spCOMPLEX and 01402 * spSEPARATED_COMPLEX_VECTORS are set, otherwise it disappears 01403 * without a trace. 01404 */ 01405 01406 void 01407 sfMultTransposed( Matrix, RHS, Solution IMAG_VECTORS ) 01408 01409 long *Matrix; 01410 RealVector Solution, RHS IMAG_VECTORS; 01411 { 01412 /* Begin `sfMultTransposed'. */ 01413 spMultTransposed( (char *)*Matrix, RHS, Solution IMAG_VECTORS ); 01414 } 01415 #endif /* MULTIPLICATION && TRANSPOSE */ 01416 01417 01418 01419 01420 01421 01422 #if DETERMINANT 01423 01424 /* 01425 * CALCULATE DETERMINANT 01426 * 01427 * This routine in capable of calculating the determinant of the 01428 * matrix once the LU factorization has been performed. Hence, only 01429 * use this routine after spFactor() and before spClear(). 01430 * The determinant equals the product of all the diagonal elements of 01431 * the lower triangular matrix L, except that this product may need 01432 * negating. Whether the product or the negative product equals the 01433 * determinant is determined by the number of row and column 01434 * interchanges performed. Note that the determinants of matrices can 01435 * be very large or very small. On large matrices, the determinant 01436 * can be far larger or smaller than can be represented by a floating 01437 * point number. For this reason the determinant is scaled to a 01438 * reasonable value and the logarithm of the scale factor is returned. 01439 * 01440 * >>> Arguments: 01441 * Matrix <input> (long *) [INTEGER] 01442 * A pointer to the matrix for which the determinant is desired. 01443 * pExponent <output> (int *) [INTEGER or INTEGER*2] 01444 * The logarithm base 10 of the scale factor for the determinant. To 01445 * find 01446 * the actual determinant, Exponent should be added to the exponent of 01447 * DeterminantReal. 01448 * pDeterminant <output> (RealNumber *) [REAL or DOUBLE PRECISION] 01449 * The real portion of the determinant. This number is scaled to be 01450 * greater than or equal to 1.0 and less than 10.0. 01451 * piDeterminant <output> (RealNumber *) [REAL or DOUBLE PRECISION] 01452 * The imaginary portion of the determinant. When the matrix is real 01453 * this pointer need not be supplied, nothing will be returned. This 01454 * number is scaled to be greater than or equal to 1.0 and less than 10.0. 01455 */ 01456 01457 #if spCOMPLEX 01458 01459 void 01460 sfDeterminant( Matrix, pExponent, pDeterminant, piDeterminant ) 01461 01462 long *Matrix; 01463 RealNumber *pDeterminant, *piDeterminant; 01464 int *pExponent; 01465 { 01466 /* Begin `sfDeterminant'. */ 01467 spDeterminant( (char *)*Matrix, pExponent, pDeterminant, piDeterminant ); 01468 } 01469 01470 #else /* spCOMPLEX */ 01471 01472 void 01473 sfDeterminant( Matrix, pExponent, pDeterminant ) 01474 01475 long *Matrix; 01476 RealNumber *pDeterminant; 01477 int *pExponent; 01478 { 01479 /* Begin `sfDeterminant'. */ 01480 spDeterminant( (char *)*Matrix, pExponent, pDeterminant ); 01481 } 01482 #endif /* spCOMPLEX */ 01483 #endif /* DETERMINANT */ 01484 01485 01486 01487 01488 01489 01490 /* 01491 * RETURN MATRIX ERROR STATUS 01492 * 01493 * This function is used to determine the error status of the given matrix. 01494 * 01495 * >>> Returned: [INTEGER or INTEGER*2] 01496 * The error status of the given matrix. 01497 * 01498 * >>> Arguments: 01499 * Matrix <input> (long *) [INTEGER] 01500 * The matrix for which the error status is desired. 01501 */ 01502 01503 int 01504 sfError( Matrix ) 01505 01506 long *Matrix; 01507 { 01508 /* Begin `sfError'. */ 01509 return spError( (char *)*Matrix ); 01510 } 01511 01512 01513 01514 01515 01516 01517 /* 01518 * WHERE IS MATRIX SINGULAR 01519 * 01520 * This function returns the row and column number where the matrix was 01521 * detected as singular or where a zero was detected on the diagonal. 01522 * 01523 * >>> Arguments: 01524 * Matrix <input> (long *) [INTEGER] 01525 * The matrix for which the error status is desired. 01526 * pRow <output> (int *) [INTEGER or INTEGER*2] 01527 * The row number. 01528 * pCol <output> (int *) [INTEGER or INTEGER*2] 01529 * The column number. 01530 */ 01531 01532 void 01533 sfWhereSingular( Matrix, Row, Col ) 01534 01535 long *Matrix; 01536 int *Row, *Col; 01537 { 01538 /* Begin `sfWhereSingular'. */ 01539 spWhereSingular( (char *)*Matrix, Row, Col ); 01540 } 01541 01542 01543 01544 01545 01546 /* 01547 * MATRIX SIZE 01548 * 01549 * Returns the size of the matrix. Either the internal or external size of 01550 * the matrix is returned. 01551 * 01552 * >>> Arguments: 01553 * Matrix <input> (long *) [INTEGER] 01554 * Pointer to matrix. 01555 * External <input> (BOOLEAN) [LOGICAL] 01556 * If External is set true, the external size , i.e., the value of the 01557 * largest external row or column number encountered is returned. 01558 * Otherwise the true size of the matrix is returned. These two sizes 01559 * may differ if the TRANSLATE option is set true. 01560 */ 01561 01562 int 01563 sfGetSize( Matrix, External ) 01564 01565 long *Matrix, *External; 01566 { 01567 /* Begin `sfGetSize'. */ 01568 return spGetSize( (char *)*Matrix, (BOOLEAN)*External ); 01569 } 01570 01571 01572 01573 01574 01575 01576 01577 01578 /* 01579 * SET MATRIX COMPLEX OR REAL 01580 * 01581 * Forces matrix to be either real or complex. 01582 * 01583 * >>> Arguments: 01584 * Matrix <input> (long *) [INTEGER] 01585 * Pointer to matrix. 01586 */ 01587 01588 void 01589 sfSetReal( Matrix ) 01590 01591 long *Matrix; 01592 { 01593 /* Begin `sfSetReal'. */ 01594 spSetReal( (char *)*Matrix ); 01595 } 01596 01597 01598 void 01599 sfSetComplex( Matrix ) 01600 01601 long *Matrix; 01602 { 01603 /* Begin `sfSetComplex'. */ 01604 spSetComplex( (char *)*Matrix ); 01605 } 01606 01607 01608 01609 01610 01611 01612 01613 01614 01615 /* 01616 * ELEMENT OR FILL-IN COUNT 01617 * 01618 * Two functions used to return simple statistics. Either the number 01619 * of total elements, or the number of fill-ins can be returned. 01620 * 01621 * >>> Arguments: 01622 * Matrix <input> (long *) [INTEGER] 01623 * Pointer to matrix. 01624 */ 01625 01626 int 01627 sfFillinCount( Matrix ) 01628 01629 long *Matrix; 01630 { 01631 /* Begin `sfFillinCount'. */ 01632 return spFillinCount( (char *)*Matrix ); 01633 } 01634 01635 01636 int sfElementCount(long *Matrix) 01637 { 01638 return spElementCount( (char *)*Matrix ); 01639 } 01640 01641 #if TRANSLATE && DELETE 01642 /* 01643 * DELETE A ROW && COLUMN FROM THE MATRIX 01644 * 01645 * Deletes a row and a column from a matrix. 01646 * 01647 * Sparse will abort if an attempt is made to delete a row or column that 01648 * doesn't exist. 01649 * 01650 * >>> Arguments: 01651 * Matrix <input> (long *) [INTEGER] 01652 * Pointer to the matrix in which the row and column are to be deleted. 01653 * Row <input> (int) [INTEGER or INTEGER*2] 01654 * Row to be deleted. 01655 * Col <input> (int) [INTEGER or INTEGER*2] 01656 * Column to be deleted. 01657 */ 01658 01659 void sfDeleteRowAndCol(long *Matrix, int *Row, int *Col) 01660 { 01661 spDeleteRowAndCol( (char *)*Matrix, *Row, *Col ); 01662 } 01663 #endif 01664 01665 #if PSEUDOCONDITION 01666 /* 01667 * CALCULATE PSEUDOCONDITION 01668 * 01669 * Computes the magnitude of the ratio of the largest to the smallest 01670 * pivots. This quantity is an indicator of ill-conditioning in the 01671 * matrix. If this ratio is large, and if the matrix is scaled such 01672 * that uncertainties in the RHS and the matrix entries are 01673 * equilibrated, then the matrix is ill-conditioned. However, a small 01674 * ratio does not necessarily imply that the matrix is 01675 * well-conditioned. This routine must only be used after a matrix 01676 * has been factored by sfOrderAndFactor() or sfFactor() and before it 01677 * is cleared by sfClear() or spInitialize(). The pseudocondition is faster 01678 * to compute than the condition number calculated by sfCondition(), but 01679 * is not as informative. 01680 * 01681 * >>> Returns: [REAL or DOUBLE PRECISION] 01682 * The magnitude of the ratio of the largest to smallest pivot used during 01683 * previous factorization. If the matrix was singular, zero is returned. 01684 * 01685 * >>> Arguments: 01686 * Matrix <input> (long *) 01687 * Pointer to the matrix. 01688 */ 01689 01690 double sfPseudoCondition(long *Matrix) 01691 { 01692 return spPseudoCondition( (char *)Matrix ); 01693 } 01694 #endif 01695 01696 #if CONDITION 01697 /* 01698 * ESTIMATE CONDITION NUMBER 01699 * 01700 * Computes an estimate of the condition number using a variation on 01701 * the LINPACK condition number estimation algorithm. This quantity is 01702 * an indicator of ill-conditioning in the matrix. To avoid problems 01703 * with overflow, the reciprocal of the condition number is returned. 01704 * If this number is small, and if the matrix is scaled such that 01705 * uncertainties in the RHS and the matrix entries are equilibrated, 01706 * then the matrix is ill-conditioned. If the this number is near 01707 * one, the matrix is well conditioned. This routine must only be 01708 * used after a matrix has been factored by sfOrderAndFactor() or 01709 * sfFactor() and before it is cleared by sfClear() or spInitialize(). 01710 * 01711 * Unlike the LINPACK condition number estimator, this routines 01712 * returns the L infinity condition number. This is an artifact of 01713 * Sparse placing ones on the diagonal of the upper triangular matrix 01714 * rather than the lower. This difference should be of no importance. 01715 * 01716 * References: 01717 * A.K. Cline, C.B. Moler, G.W. Stewart, J.H. Wilkinson. An estimate 01718 * for the condition number of a matrix. SIAM Journal on Numerical 01719 * Analysis. Vol. 16, No. 2, pages 368-375, April 1979. 01720 * 01721 * J.J. Dongarra, C.B. Moler, J.R. Bunch, G.W. Stewart. LINPACK 01722 * User's Guide. SIAM, 1979. 01723 * 01724 * Roger G. Grimes, John G. Lewis. Condition number estimation for 01725 * sparse matrices. SIAM Journal on Scientific and Statistical 01726 * Computing. Vol. 2, No. 4, pages 384-388, December 1981. 01727 * 01728 * Dianne Prost O'Leary. Estimating matrix condition numbers. SIAM 01729 * Journal on Scientific and Statistical Computing. Vol. 1, No. 2, 01730 * pages 205-209, June 1980. 01731 * 01732 * >>> Returns: [REAL or DOUBLE PRECISION] 01733 * The reciprocal of the condition number. If the matrix was singular, 01734 * zero is returned. 01735 * 01736 * >>> Arguments: 01737 * eMatrix <input> (long *) 01738 * Pointer to the matrix. 01739 * NormOfMatrix <input> (RealNumber *) [REAL or DOUBLE PRECISION] 01740 * The L-infinity norm of the unfactored matrix as computed by 01741 * spNorm(). 01742 * pError <output> (int *) [INTEGER or INTEGER*2] 01743 * Used to return error code. 01744 * 01745 * >>> Possible errors: 01746 * spSINGULAR 01747 * spNO_MEMORY 01748 */ 01749 01750 double sfCondition(long *Matrix, double NormOfMatrix, int *pError ) 01751 { 01752 return spCondition( (char *)*Matrix, *NormOfMatrix, pError ); 01753 } 01754 /* 01755 * L-INFINITY MATRIX NORM 01756 * 01757 * Computes the L-infinity norm of an unfactored matrix. It is a fatal 01758 * error to pass this routine a factored matrix. 01759 * 01760 * One difficulty is that the rows may not be linked. 01761 * 01762 * >>> Returns: [REAL or DOUBLE PRECISION] 01763 * The largest absolute row sum of matrix. 01764 * 01765 * >>> Arguments: 01766 * Matrix <input> (long *) 01767 * Pointer to the matrix. 01768 */ 01769 01770 double sfNorm(long *Matrix) 01771 { 01772 return spNorm( (char *)*Matrix ); 01773 } 01774 #endif /* CONDITION */ 01775 01776 01777 01778 01779 01780 #if STABILITY 01781 01782 /* 01783 * STABILITY OF FACTORIZATION 01784 * 01785 * The following routines are used to gauge the stability of a 01786 * factorization. If the factorization is determined to be too unstable, 01787 * then the matrix should be reordered. The routines compute quantities 01788 * that are needed in the computation of a bound on the error attributed 01789 * to any one element in the matrix during the factorization. In other 01790 * words, there is a matrix E = [e_ij] of error terms such that A+E = LU. 01791 * This routine finds a bound on |e_ij|. Erisman & Reid [1] showed that 01792 * |e_ij| < 3.01 u rho m_ij, where u is the machine rounding unit, 01793 * rho = max a_ij where the max is taken over every row i, column j, and 01794 * step k, and m_ij is the number of multiplications required in the 01795 * computation of l_ij if i > j or u_ij otherwise. Barlow [2] showed that 01796 * rho < max_i || l_i ||_p max_j || u_j ||_q where 1/p + 1/q = 1. 01797 * 01798 * The first routine finds the magnitude on the largest element in the 01799 * matrix. If the matrix has not yet been factored, the largest 01800 * element is found by direct search. If the matrix is factored, a 01801 * bound on the largest element in any of the reduced submatrices is 01802 * computed using Barlow with p = oo and q = 1. The ratio of these 01803 * two numbers is the growth, which can be used to determine if the 01804 * pivoting order is adequate. A large growth implies that 01805 * considerable error has been made in the factorization and that it 01806 * is probably a good idea to reorder the matrix. If a large growth 01807 * in encountered after using spFactor(), reconstruct the matrix and 01808 * refactor using spOrderAndFactor(). If a large growth is 01809 * encountered after using spOrderAndFactor(), refactor using 01810 * spOrderAndFactor() with the pivot threshold increased, say to 0.1. 01811 * 01812 * Using only the size of the matrix as an upper bound on m_ij and 01813 * Barlow's bound, the user can estimate the size of the matrix error 01814 * terms e_ij using the bound of Erisman and Reid. The second routine 01815 * computes a tighter bound (with more work) based on work by Gear 01816 * [3], |e_ij| < 1.01 u rho (t c^3 + (1 + t)c^2) where t is the 01817 * threshold and c is the maximum number of off-diagonal elements in 01818 * any row of L. The expensive part of computing this bound is 01819 * determining the maximum number of off-diagonals in L, which changes 01820 * only when the order of the matrix changes. This number is computed 01821 * and saved, and only recomputed if the matrix is reordered. 01822 * 01823 * [1] A. M. Erisman, J. K. Reid. Monitoring the stability of the 01824 * triangular factorization of a sparse matrix. Numerische 01825 * Mathematik. Vol. 22, No. 3, 1974, pp 183-186. 01826 * 01827 * [2] J. L. Barlow. A note on monitoring the stability of triangular 01828 * decomposition of sparse matrices. "SIAM Journal of Scientific 01829 * and Statistical Computing." Vol. 7, No. 1, January 1986, pp 166-168. 01830 * 01831 * [3] I. S. Duff, A. M. Erisman, J. K. Reid. "Direct Methods for Sparse 01832 * Matrices." Oxford 1986. pp 99. 01833 */ 01834 01835 /* 01836 * LARGEST ELEMENT IN MATRIX 01837 * 01838 * >>> Returns: [REAL or DOUBLE PRECISION] 01839 * If matrix is not factored, returns the magnitude of the largest element in 01840 * the matrix. If the matrix is factored, a bound on the magnitude of the 01841 * largest element in any of the reduced submatrices is returned. 01842 * 01843 * >>> Arguments: 01844 * Matrix <input> (long *) [INTEGER] 01845 * Pointer to the matrix. 01846 */ 01847 01848 double sfLargestElement(long *Matrix) 01849 { 01850 return spLargestElement( (char *)Matrix ); 01851 } 01852 /* 01853 * MATRIX ROUNDOFF ERROR 01854 * 01855 * >>> Returns: [REAL or DOUBLE PRECISION] 01856 * Returns a bound on the magnitude of the largest element in E = A - LU. 01857 * 01858 * >>> Arguments: 01859 * Matrix <input> (long *) [INTEGER] 01860 * Pointer to the matrix. 01861 * Rho <input> (RealNumber *) [REAL or DOUBLE PRECISION] 01862 * The bound on the magnitude of the largest element in any of the 01863 * reduced submatrices. This is the number computed by the function 01864 * spLargestElement() when given a factored matrix. If this number is 01865 * negative, the bound will be computed automatically. 01866 */ 01867 01868 double sfRoundoff(long *Matrix, double Rho) 01869 { 01870 return spRoundoff( (char *)*Matrix, *Rho ); 01871 } 01872 #endif 01873 #endif /* FORTRAN */