/*Translated by FOR_C, v3.4.2 (-), on 07/09/115 at 08:31:22 */ /*FOR_C Options SET: ftn=u io=c no=p op=aimnv s=dbov str=l x=f - prototypes */ #include #include "fcrt.h" #include "sspge.h" #include #include /* PARAMETER translations */ #define IPBDIM 2 #define IPLAST 5 #define LTXTAC 54 #define LTXTAD 137 #define LTXTAE 182 #define LTXTAF 225 #define LTXTAG 282 #define LTXTAH 316 #define LTXTAI 370 #define LTXTAJ 413 #define MECONT 50 #define MEEMES 52 #define MENTXT 23 #define MERET 51 #define METEXT 53 /* end of PARAMETER translations */ void /*FUNCTION*/ sspge( long n, long ispec[], long ia[], float a[], float b[], float opt[]) { LOGICAL32 cond, det, trans; long int i, ib, icp, idat[4], il, incb, ip, ipivot, ir, j, k, kb, kk, l, lastr, lc, lcfw, lclisc, lcoli, lcolsz, lcp, lfree, lnext, lr, lrowl, lrowli, lrowsz, lrswsc, m, mrowsz, nb, ni, np; float fdat[2], pval, tp, tp1, tp2, tp3; void *_p0; static char mtxtaa[3][210]={"SSPGE$BN = $I, IA dim. = $I, A dim. = $I, Options: $EWith $I columns yet to factor, our best pivot is $F$NThe matrix appears singular.$EColumn $I is all 0, the matrix is singular.$ERow$ $I is all 0, the matrix", " is singular.$EWith $I columns yet to factor we have run out of space.$EThe dimension provided is .le. 0$EThe last option index you provided is not supported.$EThe column$ with index $I has $I entries?$EIn col", "umn $I, row index $I follows row index $I?$NRow indexes must be increasing and in the inteval [1,N].$EA bug in SSPGE -- Could not find unused column.$EA bug$ in SSPGE -- Pivot row not found where expected.$E "}; static char mtxtab[1][14]={"$I $I $I; $B"}; static char mtxtac[1][8]={"$I.$N$B"}; static long mloc[8]={LTXTAC,LTXTAD,LTXTAE,LTXTAF,LTXTAG,LTXTAH, LTXTAI,LTXTAJ}; static long macter[5]={MEEMES,0,0,0,MECONT}; static long mactop[2]={METEXT,MECONT}; static long macttx[4]={MENTXT,0,METEXT,MERET}; /* OFFSET Vectors w/subscript range: 1 to dimension */ float *const A = &a[0] - 1; float *const B = &b[0] - 1; float *const Fdat = &fdat[0] - 1; long *const Ia = &ia[0] - 1; long *const Idat = &idat[0] - 1; long *const Ispec = &ispec[0] - 1; long *const Macter = &macter[0] - 1; long *const Mactop = &mactop[0] - 1; long *const Macttx = &macttx[0] - 1; long *const Mloc = &mloc[0] - 1; float *const Opt = &opt[0] - 1; /* end of OFFSET VECTORS */ /* Time-stamp: <2015-07-08 15:34:13 m> * Copyright (c) 2006, Math a la Carte, Inc. */ /* Solve the system of linear equations, A X = B, return X in B. *>> 2015-07-08 SSPGE Krogh -- Probably this code is worthless *>> 2006-03-12 SSPGE Krogh Initial code */ /*--S replaces "?": ?SPGE, ?MESS */ /* ************************* Arguments ********************************** */ /* A The input matrix and associated working store. * IA Array for holding data defining the sparity structure and also * used for working storage. IA(1) contains the number of nonzeros in * column 1, the next IA locations contain the row indexes for this * column in increasing order. This is followed by the number of * nonzeros in column 2, etc. Additional working space is needed to * hold data as we factor the matrix. The data in A is stored at the * corresponding locations in A, with anything being stored in A where * the column counts are stored in IA. * ISPEC Used as follows. * ISPEC(1) Must be 0 if the matrix is not factored and the value * returned on the first call must be preserved on later calls when * the matrix already factored. If the returned value is < 0, * then the matrix was not factored and the value is the negative of * the values mentioned in Error message indexes below. * ISPEC(2) Declared dimension of IA() (or at least all that you want * this routine to know about). Typically this should be the same as * ISPEC(3) + 4N. * ISPEC(3) [in] Declared dimension of A(). This number should be the * space required for the original matrix + the maximal space expected * for the same type of representation of the factored matrix + 2N. * ISPEC(4) The lesser of the amount of space in A or IA that was not * needed. * ISPEC(5) Start of options, * OPT Used for return of optional values. * OPT(1) The reciprocal of the condition number if requested. * OPT(2) If requested the determinant is OPT(2) * 10 ** OPT(3) * OPT(3) The exponent for the determinant. (0 if the determinant is * not to big or too small. * Error message indexes. * 1 Matrix appears singular. * 2 Matrix has an empty column. * 3 Matrix has an empty row. * 4 No more space. * 5 N .le. 0 * 6 An unknown option * 7 Problem with column size * 8 Problem with row indexes * 9 Bug -- Can't find an unused column * 10 Bug -- Pivot row not found where expected * * ************************ Other Variables ***************************** * * COND Usually .false., set =.true. if condition estimate wanted. * DET Usually .false., set=.true. if determinant is wanted. * I Temporary index * IB Base address for the current column in B. * IC Temporary index of a column. * ICP Initial index of current pivot column. * IL Temporary index. * INCB Increment between columns of B. * IP Physical row index. * IPIVOT Logical index of the pivot column. * IR Temporary index of a row. * J Temporary index * K Temporary index. * KK Temporary index. * KB Index of best choice for a pivot and counter on columns in B. * L Temporary index. * LASTR Index of last logical row used for selecing the next column. * LC Location of a column. * LCFW Base address in IA mapping logical index of a column to its * factored location. Corresponding space in A is used to accumulate * transformations into the next pivot column. * LCLISC Base address in IA mapping logical index of a column to its * initial index. Values here are set negative if this column has * been used to adjust row counts, and is set to 0 if the column * is used as a pivot column. Also used as the base address for * the column scaling. * LCOLI Base address in IA mapping intial index of a column to its * location. * LCP Location of original column data for the pivot column. * LCOLSZ Base address for adjusted size of coluumns. These are stored * in initial column order. * LFREE Address of the last available free space in IA. * LNEXT Base location where next factored data is saved. * LR Location of a row. * LROWL Base address for pointers to the start of list of initial * columns for each row. These are stored in initial row order. * The pointers point to the number of nonzero columns in the row * followed by the list of column entries. * LROWLI Base address for mapping logical row indexes to initial * indexes. * LRSWSC Base address in IA for tracking row swaps. Also used for the * base address of the row scaling. The scales are stored in iniitial * row order. * LROWSZ Base address for adjusted size of rows. These are stored * in logical row order. */ /* MROWSZ The smallest adjusted row size. * NB Number of columns in B. * NI Internal value of N. * NP Index of a pivot in U. * PVAL Value in pivot location (if no swap) when factoring. Also * used for reciprocal of the L1 condition number of the orig. matrix. * TP,TP1,TP2,TP3 Temporary variables.. * TRAN Usuall .false., set = .true. if solving U^T x = b. */ /* Notes on permutation: * When permuting the mapping from logical to iniital indexes, we need * to know the two logical indexes and then we swap the contents of * IA at these two logical indexes. * When permuting the mapping from initial to logical indexes, we need * to know the two intial indexes and then we swap the contents of * IA at these tow initial indexes. * It is frequently the case that we need one mapping to get the * indexes for permuting the other. * Given the indices to swap, say IL1, IL2, and LI1, and LI2, where * the IL indexes are initial indexes which we map to logical indexes, * and the LI indexes are logical indexes whichg map to initial indexes. * Then New IL1 mapping is LI2, New IL2 mapping is LI1, and * New LI1 mapping is IL2, New LI2 mapping is IL1. * Once IL1, IL2, LI1, LI2 have been determined new assignments are * [IL1] = LI2 [IL2] = LI1 [LI1] = IL2 [LI2= IL1 * where [IL1] means location IL1 in the permutation mapping initial * indexes to logical, etc. * * Notes on storage in IA and A. * First Location Used * Contents (IA / A ) IA A * Input matrix 1 1 * Column: Logical => Initial Perm/Scaling LCLISC+1 LCLISC+1 * Row swap history/Row Scaling LRSWSC+1 LRSWSC+1 * Map logical col => factored loc. / Various LCFW+1 LCFW+1 * * Place to store rest of columns of factored part LNEXT+1 LNEXT+1 * End of space available for factored part. LFREE LFREE * Data below only uses space while factoring the matrix. * Map initial row index to data for that row. LROWL+1 * Adjusted row sizes LROWSZ+1 * Initial size of columns LCOLSZ+1 * Pointer to columns of the original matrix LCOLI+1 * Map logical to initial row indexes LROWLI+1 * */ /* A factored column has the form * s The space used by the column. (location LNEXT+1 or IA(LCFW+J)) * p The location of the pivot is p. (location LNEXT+2 or IA(LCFW+J+1)) * p - 2 locations containing a column of U. (locations LNEXT+3 to p) * The reciprocal of the pivot. (location p) * s - (p-LNEXT) - 1 locations for the subdiagonal entries of the * L factors (locations p+1 to LNEXT+s-1) * Except for the first two entries, IA contains the logical row indexes * for U and L, and A contains the corresponding values. * * ************************** Declarations ****************************** */ /* ********* Error message text *************** *[Last 2 letters of Param. name] [Text generating message.] *AA SSPGE$B *AB N = $I, IA dim. = $I, A dim. = $I, Options: $E *AC With $I columns yet to factor, our best pivot is $F$N * The matrix appears singular.$E *AD Column $I is all 0, the matrix is singular.$E *AE Row $I is all 0, the matrix is singular.$E *AF With $I columns yet to factor we have run out of space.$E *AG The dimension provided is .le. 0$E *AH The last option index you provided is not supported.$E *AI The column with index $I has $I entries?$E *AJ In column $I, row index $I follows row index $I?$N * Row indexes must be increasing and in the inteval [1,N].$E *AK A bug in SSPGE -- Could not find unused column.$E *AL A bug in SSPGE -- Pivot row not found where expected.$E * $ *AM $I $I $I; $B * $ *AN $I.$N$B */ /* *********** End of generated error message data ********** */ /*++ Default DEBUG = 0 */ /*++ Code for DEBUG > 0 is inactive * integer MEFSPV, MEFVEC, MEIVEC, METDIG, MEMDA1, MEMDA2 * parameter (MEFSPV=54, MEFVEC=61, MEIVEC=57, METDIG=22, * 1 MEMDA1=27, MEMDA2=28) * integer MACTSP(10), MACTFV(6), MACTSV(6), MACTIV(4), MACTXX(2) * character SPCOL(1)*32, SPROW(1)*32, BVEC(1)*4, COLSCA(1)*16, * 1 ROWSCA(1)*16, PCOLLI(1)*36, ROWSZ(1)*12, COLSZ(1)*12, * 2 COLLOC(1)*11, PROB(1)*46, COLU(1)*60, RMROW(1)*56 * integer DIGS * parameter (DIGS=8) *c * data SPCOL(1) / ' Col $M, $M Row Entries$B' / * data SPROW(1) / ' Row $M, $M Col Entries$B' / * data BVEC(1) / 'B:$B'/ * data COLSCA(1) / 'Column Scales:$B'/ * data ROWSCA(1) / ' Row Scales:$B'/ * data PCOLLI(1) / 'Logical => Initial Column mapping:$E' / * data ROWSZ(1) / 'Row sizes:$E' / * data COLSZ(1) / 'Col sizes:$E' / * data COLLOC(1) /'Col locs:$E' / * data RMROW(1) / * 1 'Col $I removed from row $I, (logical row $I), $I left.$E' / * data PROB(1) / 'N = $I, Nonzeros = $I, Fraction nonzero = $F$E'/ * data COLU(1) / * 1 'Col. $I factored to col. $I with pivot $F in init. row $I.$E' / *c 1 2 3 4 5 6 7 8 9 10 * data MACTSP/MEMDA1,1,MEMDA2,1,METEXT,METDIG,DIGS,MEFSPV,1,MERET/ *c 1 2 3 4 5 6 * data MACTFV / METEXT, METDIG, DIGS, MEFVEC, 0, MERET / * data MACTSV / METEXT, METDIG, 4, MEFVEC, 0, MERET / * data MACTIV / METEXT, MEIVEC, 0, MERET / * data MACTXX / METEXT, MERET / * MACTIV(3) = N * MACTFV(5) = N * MACTSV(5) = N *++ End */ /* integer DIGS * parameter (DIGS=8) * integer MEFSPV, MEFVEC, MEIVEC, METDIG, MEMDA1, MEMDA2 * parameter (MEFSPV=54, MEFVEC=61, MEIVEC=57, METDIG=22, * 1 MEMDA1=27, MEMDA2=28) * character NEWROW(1)*27 * integer MACTNR(10) *c 1 2 3 4 5 6 7 8 9 10 * data MACTNR/MEMDA1,1,MEMDA2,1,METEXT,METDIG,DIGS,MEFSPV,1,MERET/ * data NEWROW(1) / 'Row at $M, $M Row Entries$B' / */ ni = n; if (ni <= 0) { Macter[3] = 5; goto L_2000; } nb = 1; incb = 1; trans = FALSE; cond = FALSE; det = FALSE; /* Options checked here */ k = 4; L_10: ; k += 1; j = Ispec[k]; if (j > 1) { if (j == 2) { /* Setting the number of columns in B. */ nb = Ispec[k + 1]; incb = Ispec[k + 2]; k += 2; } else if (j == 3) { cond = TRUE; } else if (j == 4) { det = TRUE; } else { /* Unknnown option */ Macter[3] = 6; goto L_2000; } goto L_10; } trans = j == 1; if (Ispec[1] == 0) { lc = 1; l = 0; for (i = 1; i <= n; i++) { /* Get size of user data. */ m = Ia[lc + l]; if ((m <= 0) || (m > ni)) { /* Bad column size */ Macter[3] = 7; Idat[1] = i; Idat[2] = m; goto L_2000; } if (l == 0) { for (k = lc + 1; k <= (lc + m); k++) { if (A[k] == 0.e0) goto L_80; } lc += m + 1; goto L_110; } /* Matrix has 0 stored, remove them. */ L_80: ; kb = lc + l; kk = m; for (k = kb + 1; k <= (kb + m); k++) { if (A[k] != 0) { A[k - l] = A[k]; Ia[k - l] = Ia[k]; } else { l += 1; kk -= 1; } } Ia[lc] = kk; lc += kk + 1; L_110: ; } } else { lc = Ispec[1]; } /* Set pointers for internal vectors. */ lclisc = lc - 1; lrswsc = lc + ni; lcfw = lrswsc + ni; lnext = lcfw + ni; /* Go directly to solve if matrix is factored. */ if (Ispec[1] > 0) goto L_750; Ispec[1] = lc; lrowli = Ispec[2] - ni; lcoli = lrowli - ni; lcolsz = lcoli - ni; lrowsz = lcolsz - ni; lrowl = lrowsz - ni; Ia[lrowl + 1] = min( lrowl, Ispec[3] ) - lc; Ia[Ia[lrowl + 1]] = 0; if (lnext + ni >= Ia[lrowl + 1]) { /* Insufficient memory */ Macter[3] = 4; Idat[1] = ni; goto L_2000; } /* Initiialize scaling (Later an option to avoid this may be added.) */ for (i = lc; i <= lcfw; i++) { A[i] = 0.e0; } for (i = lrowsz + 1; i <= lcolsz; i++) { Ia[i] = 0; } /* Get sizes, check indexes. */ lc = 1; for (j = 1; j <= ni; j++) { Ia[lrowli + j] = j; Ia[lclisc + j] = j; Ia[lcoli + j] = lc; A[lcfw + j] = 0.e0; i = 0; for (k = lc + 1; k <= (lc + Ia[lc]); k++) { if ((Ia[k] <= i) || (Ia[k] > ni)) { /* Bad row indexes */ Macter[3] = 8; Idat[1] = j; Idat[2] = Ia[k]; Idat[3] = i; goto L_2000; } i = Ia[k]; Ia[lrowsz + i] += 1; /* Column scaling */ A[lclisc + j] = fmaxf( A[lclisc + j], fabsf( A[k] ) ); } lc = k; } /* Get row data locations and index of shortest row. */ mrowsz = Ia[lrowsz + 1]; lastr = 1; for (i = 2; i <= ni; i++) { Ia[lrowl + i] = Ia[lrowl + i - 1] + Ia[lrowsz + i - 1] + 1; if (Ia[lrowsz + i] < mrowsz) { mrowsz = Ia[lrowsz + i]; lastr = i; } Ia[Ia[lrowl + i]] = 0; } /* Fix scaling and get data indexes ordered by row then column. */ for (j = 1; j <= ni; j++) { if (A[lclisc + j] == 0.e0) { /* Column is all 0. */ Macter[3] = 2; Idat[1] = j; Macter[2] = 28; goto L_2010; } A[lclisc + j] = 1.e0/A[lclisc + j]; l = Ia[lcoli + j]; Ia[lcolsz + j] = Ia[l]; for (k = l + 1; k <= (l + Ia[l]); k++) { i = Ia[k]; il = Ia[lrowl + i]; Ia[il] += 1; Ia[il + Ia[il]] = j; A[il + Ia[il]] = A[k]*A[lclisc + j]; /* Update row scale. */ A[lrswsc + i] = fmaxf( A[lrswsc + i], fabsf( A[il + Ia[il]] ) ); } } lfree = Ia[lrowl + 1] - ni; /*++ Code for DEBUG > 0 is inactive * IDAT(1) = NI * IDAT(2) = LCLISC - NI * FDAT(1) = real(IDAT(2)) / real(NI*NI) * call SMESS(MACTXX, PROB, IDAT, FDAT) *++ End *++ Code for DEBUG > 9 is inactive *c Print the initial matrix * K = 1 * do 380 I = 1, N * MACTSP(2) = I * MACTSP(4) = IA(K) * MACTSP(9) = IA(K) * call SMESS(MACTSP, SPCOL, IA(K+1), A(K+1)) * K = K + IA(K) + 1 * 380 continue *++ End *++ Code for DEBUG > 7 is inactive *c Print the initial B * call SMESS(MACTFV, BVEC, IDAT, B) *++ End *++ Code for DEBUG > 4 is inactive *c Print the scaling * call SMESS(MACTSV,COLSCA, IDAT, A(LCLISC+1)) * call SMESS(MACTSV,ROWSCA, IDAT, A(LRSWSC+1)) *++ End *++ Code for DEBUG > 4 is inactive *c Print permutations and row sizes and column locations * call MESS(MACTIV,PCOLLI, IA(LCLISC+1)) * call MESS(MACTIV,ROWSZ, IA(LROWSZ+1)) * call MESS(MACTIV,COLLOC, IA(LCOLI+1)) *++ End */ if (2*lcfw + ni >= lfree) { /* Insuffient memory */ Macter[3] = 4; Idat[1] = ni; goto L_2000; } for (i = 1; i <= ni; i++) { /* Save reciprocal of the row scale factors */ if (A[lrswsc + i] == 0.e0) { /* Row is all 0. */ Macter[3] = 3; Idat[1] = i; Macter[2] = 28; goto L_2010; } A[lrswsc + i] = 1.e0/A[lrswsc + i]; } /* ********************** Do the Factorization ************************** */ /*??? Need some checks various places for out of memory. */ lr = Ia[lrowl + lastr]; lastr = ni; for (ipivot = ni; ipivot >= 1; ipivot--) { /* Pick the pivot column */ L_400: ; /*++ Code for DEBUG > 1 is inactive *c Print the currently selected row * MACTSP(2) = IA(LROWLI+LASTR) * MACTSP(4) = IA(LR) * MACTSP(9) = IA(LR) * call SMESS(MACTSP, SPROW, IA(LR+1), A(LR+1)) *++ End */ tp = 1.e30; kb = 0; j = Ia[lr]; l = 1000000000; for (k = lr + 1; k <= (lr + j); k++) { L_405: ; if (Ia[lcoli + Ia[k]] == 0) { /*++ Code for DEBUG > 1 is inactive *c Print col. removed from row. * IDAT(1) = IA(K) * IDAT(2) = IA(LROWLI+LASTR) * IDAT(3) = LASTR * IDAT(4) = IA(LR) - 1 * call MESS(MACTSP, RMROW, IDAT) *++ End */ Ia[k] = Ia[lr + Ia[lr]]; A[k] = A[lr + Ia[lr]]; Ia[lr] -= 1; Ia[lrowsz + Ia[lrowli + lastr]] = Ia[lr]; if (k <= lr + Ia[lr]) goto L_405; goto L_415; } kk = Ia[lcolsz + Ia[k]]; if (lastr > ipivot) { if (l > kk) { l = kk; kb = k; } } else if (fabsf( A[k] )*tp > (float)( kk )) { tp = (float)( kk )/fabsf( A[k] ); kb = k; } if (k > lr + Ia[lr]) goto L_415; } L_415: ; if (kb == 0) { /* We need to pick another row to get the column */ lastr -= 1; if (lastr <= ipivot) { kb = 1000000000; for (i = 1; i <= ipivot; i++) { ir = Ia[lrowli + i]; k = Ia[lrowsz + ir]; if (k < kb) { if (k != 0) { lastr = i; kb = k; } } } if (kb == 1000000000) { lastr = ipivot; for (icp = 1; icp <= ni; icp++) { if (Ia[lcoli + icp] != 0) goto L_428; } /* Can't find an unused column -- a bug. */ Macter[3] = 9; goto L_2000; } } lr = Ia[lrowl + Ia[lrowli + lastr]]; lastr = max( lastr, ipivot ); goto L_400; } else { icp = Ia[kb]; } L_428: ; lcp = labs( Ia[lcoli + icp] ); /*++ Code for DEBUG > 1 is inactive *c Print the currently selected column * MACTSP(2) = ICP * MACTSP(4) = IA(LCP) * MACTSP(9) = IA(LCP) * call SMESS(MACTSP, SPCOL, IA(LCP+1), A(LCP+1)) *++ End * Save logical to initial mapping for factored columns. */ Ia[lclisc + ipivot] = icp; /* Copy in the original column data */ for (i = lcp + 1; i <= (lcp + Ia[lcp]); i++) { A[lcfw + Ia[i]] = A[i]*A[lrswsc + Ia[i]]*A[lclisc + icp]; } /* Apply the transformations *??? This loop could be stopped one sooner, and then do * a last one storing directly to LNEXT?? */ for (j = ni; j >= (ipivot + 1); j--) { lc = Ia[lcfw + j]; np = Ia[lc + 1]; /* Row swap */ tp = A[lcfw + Ia[lrswsc + j]]; A[lcfw + Ia[lrswsc + j]] = A[lcfw + j]; A[lcfw + j] = tp; if (tp != 0.e0) { /* There is a transformation to do */ tp *= -A[np]; for (i = lc + 2; i <= (np - 1); i++) { A[lcfw + Ia[i]] += tp*A[i]; } A[lcfw + j] = tp; /*++ Code for DEBUG > 8 is inactive * print '(''IPIVOT,J:'', 2I6)', IPIVOT, J * call DVECPR(A(LCFW+1), N, 'New Col:', 0, 0, 8) *++ End */ } } Ia[lcfw + ipivot] = lnext + 1; k = lnext + 2; pval = A[lcfw + ipivot]; tp = 1.e30; for (i = 1; i <= ipivot; i++) { /* Copy U data to factored space and select the new pivot. */ if (A[lcfw + i] != 0.e0) { k += 1; Ia[k] = i; A[k] = A[lcfw + i]; A[lcfw + i] = 0.e0; tp1 = powif((float)( Ia[lrowsz + Ia[lrowli + i]] ) + .125e0,2); j = Ia[lrowl + Ia[lrowli + i]]; l = Ia[j]; for (kk = j + 1; kk <= (j + l); kk++) { if (Ia[kk] == icp) { /*++ Code for DEBUG > 1 is inactive *c Print col. removed from row. * IDAT(1) = ICP * IDAT(2) = IA(LROWLI+I) * IDAT(3) = I * IDAT(4) = IA(J) - 1 * call MESS(MACTSP, RMROW, IDAT) *++ End */ Ia[kk] = Ia[j + Ia[j]]; A[kk] = A[j + Ia[j]]; Ia[j] -= 1; goto L_495; } } L_495: ; Ia[lrowsz + Ia[lrowli + i]] = Ia[j]; if (fabsf( A[k] )*tp > tp1) { tp = tp1/fabsf( A[k] ); kb = k; } } } if (tp >= 1.e18) { /* Matrix appears to be singular. */ if (tp == 1.e30) { Fdat[1] = 0.e0; } else { Fdat[1] = A[kb]; } Idat[1] = ipivot; Macter[3] = 1; Opt[1] = 0.e0; Ispec[4] = -ipivot; Macter[2] = 28; goto L_2010; } ip = Ia[kb]; /* Get pivot stored where it should be */ if (ip == ipivot) { Ia[lrswsc + ipivot] = ipivot; } else { Ia[lrswsc + ipivot] = ip; j = Ia[lrowli + ipivot]; Ia[lrowli + ipivot] = Ia[lrowli + Ia[kb]]; Ia[lrowli + Ia[kb]] = j; if (pval != 0.e0) { /* Just an exchange */ if (Ia[k] != ipivot) { /* Pivot row not found where expected -- a bug. */ Macter[3] = 10; goto L_2000; } A[k] = A[kb]; A[kb] = pval; } else { /* Old data needs to move up. */ tp = A[kb]; for (i = kb; i <= (k - 1); i++) { Ia[i] = Ia[i + 1]; A[i] = A[i + 1]; } A[i] = tp; Ia[k] = ipivot; } } if (lastr == ipivot) lr = Ia[lrowl + Ia[lrowli + ipivot]]; Ia[lnext + 2] = k; /* Save reciprocal of the pivot. */ A[k] = 1.e0/A[k]; for (i = ipivot + 1; i <= ni; i++) { /* Copy L data to factored space. */ if (A[lcfw + i] != 0.e0) { k += 1; Ia[k] = i; A[k] = A[lcfw + i]; A[lcfw + i] = 0.e0; } } Ia[lnext + 1] = k - lnext; Ia[lcoli + icp] = 0; /*++ Code for DEBUG > 1 is inactive *c Print data about the pivot. * IDAT(1) = IA(LCLISC+IPIVOT) * IDAT(2) = IPIVOT * IDAT(3) = IA(LROWLI+IPIVOT) * FDAT(1) = 1.E0 / A(IA(LNEXT+2)) * call SMESS(MACTXX, COLU,IDAT, FDAT) * call MESS(MACTIV, ROWSZ, IA(LROWSZ+1)) * call MESS(MACTIV, COLSZ, IA(LCOLSZ+1)) *++ End */ /*++ Code for DEBUG > 6 is inactive * call MESS(MACTIV,PCOLLI, IA(LCLISC+1)) * MACTSP(2) = IPIVOT * MACTSP(4) = IA(LNEXT+1) - 2 * MACTSP(9) = MACTSP(4) * call SMESS(MACTSP, SPCOL, IA(LNEXT+3), A(LNEXT+3)) * call MESS(MACTIV,COLLOC, IA(LCFW+1)) *++ End */ lnext = k; if (lnext >= lfree) { /* Insufficient space */ Macter[3] = 4; Macter[2] = 16; Idat[1] = ipivot; Ispec[4] = -ni; goto L_2010; } } /* ****************** Compute A^{-1} b and A^{-T} b ******************** * * We have A x = b => (P_r A P_c) P_c^{-1} x = P_r b * = (U L^{-1}) P_c^{-1} x = P_r b * Thus x = P_c L U^{-1} P_r b, where * P_r = permutations mapping initial row indexes to locgical row indexes * P_c = permutations mapping initial col indexes to locgical col indexes * */ /* First get b */ L_750: ; if (nb == 0) goto L_1000; ib = 0; kb = 1; if (nb < 0) { nb = ni; for (j = 1; j <= ni; j++) { for (i = ib + 1; i <= (ib + ni); i++) { B[i] = 0.e0; } B[ib + j] = 1.e0; ib += incb; } ib = 0; } if (trans) { /* Solving the transposed system A^T x = b */ L_810: ; for (i = 1; i <= ni; i++) { /* Store the B vector in logical column order and scale.. */ A[lcfw + i] = B[ib + Ia[lclisc + i]]*A[lclisc + Ia[lclisc + i]]; } /* Get L^T P_c b */ for (j = ni; j >= 1; j--) { tp = A[lcfw + j]; k = Ia[lcfw + j]; for (i = Ia[k + 1] + 1; i <= (k + Ia[k] - 1); i++) { tp += A[i]*A[lcfw + Ia[i]]; } A[lcfw + j] = tp; } /*++ Code for DEBUG > 8 is inactive * call DVECPR(A(LCFW+1), N, 'After L^{T}:', 0, 0, 8) *++ End * Now U^{-T} of the above. */ for (j = 1; j <= ni; j++) { k = Ia[lcfw + j]; tp = A[lcfw + j]; for (i = k + 2; i <= (Ia[k + 1] - 1); i++) { tp += -A[i]*A[lcfw + Ia[i]]; } tp *= A[i]; /* Row swap */ l = Ia[lrswsc + j]; A[lcfw + j] = A[lcfw + l]; A[lcfw + l] = tp; } /*++ Code for DEBUG > 8 is inactive * call DVECPR(A(LCFW+1), N, 'After U^{-T}:', 0, 0, 8) *++ End * Get stored back in B (with scaling. */ for (i = 1; i <= ni; i++) { B[ib + i] = A[lcfw + i]*A[lrswsc + i]; } if (kb < nb) { kb += 1; ib += incb; goto L_810; } } else { /* Solving A x = b */ L_910: ; for (i = 1; i <= ni; i++) { A[lcfw + i] = B[ib + i]*A[lrswsc + i]; } /* The backsolve using w = U^{-1} b */ for (j = ni; j >= 1; j--) { lc = Ia[lcfw + j]; np = Ia[lc + 1]; tp = A[lcfw + Ia[lrswsc + j]]; A[lcfw + Ia[lrswsc + j]] = A[lcfw + j]; A[lcfw + j] = tp*A[np]; for (i = lc + 2; i <= (np - 1); i++) { A[lcfw + Ia[i]] += -A[lcfw + j]*A[i]; } } /*++ Code for DEBUG > 8 is inactive * call DVECPR(A(LCFW+1), N, 'After U^{-1}:', 0, 0, 8) *++ End * Now x = L w */ for (j = 1; j <= (ni - 1); j++) { lc = Ia[lcfw + j]; for (i = Ia[lc + 1] + 1; i <= (lc + Ia[lc] - 1); i++) { A[lcfw + Ia[i]] += A[lcfw + j]*A[i]; } } /*++ Code for DEBUG > 8 is inactive * call DVECPR(A(LCFW+1), N, 'Perm. X:', 0, 0, 8) *++ End * Now get x in correct permuted order */ for (i = 1; i <= ni; i++) { B[ib + Ia[lclisc + i]] = A[lcfw + i]*A[lclisc + Ia[lclisc + i]]; } /*++ Code for DEBUG > 7 is inactive * call DVECPR(A(LNEXT+1), N, ' X:', 0, 0, 8) *++ End */ if (kb < nb) { kb += 1; ib += incb; goto L_910; } } L_1000: ; if (cond) { /* Get estimate of reciprocal of the condition number. We ignore column * permutations as column permuatations don't affect condition number. * Get the condition number of the original scaled matrix. */ pval = 0.e0; lc = 1; for (j = 1; j <= ni; j++) { tp = 0.e0; for (i = lc + 1; i <= (lc + Ia[lc]); i++) { tp += fabsf( A[i] )*A[lrswsc + Ia[i]]; } lc = i; pval = fmaxf( pval, tp*A[lclisc + i] ); } /* Solve U y = +/- 1, each choice picks biggest norm. */ for (j = lcfw + 1; j <= (lcfw + ni); j++) { A[j] = 0.e0; } for (j = ni; j >= 1; j--) { lc = Ia[lcfw + j]; np = Ia[lc + 1]; /* Row "swap" */ tp = A[lcfw + Ia[lrswsc + j]]; A[lcfw + Ia[lrswsc + j]] = A[lcfw + j]; /* Want to try both +1 and -1 here. */ tp2 = A[np]*(tp + 1.e0); tp3 = A[np]*(tp - 1.e0); tp = fabsf( tp2 ); tp1 = fabsf( tp3 ); for (i = lc + 2; i <= (np - 1); i++) { tp += fabsf( A[lcfw + Ia[i]] - tp2*A[i] ); tp1 += fabsf( A[lcfw + Ia[i]] - tp3*A[i] ); } if (tp > tp1) { tp = tp2; } else { tp = tp3; } A[lcfw + j] = tp; for (i = lc + 2; i <= (np - 1); i++) { A[lcfw + Ia[i]] += -tp*A[i]; } } /* Multiply this solution times L */ for (j = 1; j <= (ni - 1); j++) { lc = Ia[lcfw + j]; for (i = Ia[lc + 1] + 1; i <= (lc + Ia[lc] - 1); i++) { A[lcfw + Ia[i]] += A[lcfw + j]*A[i]; } } /* Now doing a backsolve with A^T. * Get max norm of previous result and multiply by L^T */ tp1 = 0.e0; for (j = ni; j >= 1; j--) { tp = A[lcfw + j]; tp1 = fmaxf( tp1, fabsf( tp ) ); k = Ia[lcfw + j]; for (i = Ia[k + 1] + 1; i <= (k + Ia[k] - 1); i++) { tp += A[i]*A[lcfw + Ia[i]]; } A[lcfw + j] = tp; } /* Now U^{-T} of the above. */ tp2 = 0.e0; for (j = 1; j <= ni; j++) { k = Ia[lcfw + j]; tp = A[lcfw + j]; for (i = k + 2; i <= (Ia[k + 1] - 1); i++) { tp += -A[i]*A[lcfw + Ia[i]]; } tp *= A[i]; l = Ia[lrswsc + j]; A[lcfw + j] = A[lcfw + l]; A[lcfw + l] = tp; tp2 = fmaxf( tp2, fabsf( A[lcfw + j] ) ); } Opt[1] = tp1/(pval*tp2); } if (det) { /* Compute the determinate */ Opt[3] = 0.e0; /* First get sign change due to permutations. */ tp = 1; for (i = 1; i <= ni; i++) { /* First for column swaps. */ if (Ia[lclisc + i] < 0) { Ia[lclisc + i] = -Ia[lclisc + i]; } else if (Ia[lclisc + i] != i) { j = i; L_1400: ; Ia[lclisc + j] = -Ia[lclisc + j]; j = -Ia[lclisc + j]; if (j != i) { tp = -tp; goto L_1400; } } /* Row swaps are easier */ if (Ia[lrswsc + i] != i) tp = -tp; } for (j = 1; j <= ni; j++) { k = Ia[Ia[lcfw + j] + 1]; tp *= A[k]*(A[lrswsc + j]*A[lclisc + j]); if (fabsf( tp ) > 1.e25) { tp *= 1.e-25; Opt[3] -= 25.e0; } else if (fabsf( tp ) < 1.e-25) { tp *= 1.e25; Opt[3] += 25.e0; } } Opt[2] = 1.e0/tp; } Ispec[4] = lfree - lnext - 1; return; /* *********************** Error processing ***************************** * */ L_2000: ; Macter[2] = 88; L_2010: ; Macttx[2] = Mloc[Macter[3]]; /* Output start of error message */ Ispec[1] = ni; mess( macter, (char*)mtxtaa,210, &Ispec[1] ); /* Output of options */ k = 5; L_2020: ; j = 7; i = 1; if (Ispec[k] == IPBDIM) { j = 1; i = 3; } if ((Ispec[k] > 3) && (Ispec[k] <= IPLAST)) { mess( mactop, (_p0=ntstr(mtxtab[0]+j - 1,-j + 14)),-j + 15, &Ispec[k] ); f_subscpy( mtxtab[0], j - 1, -1, 13, _p0 ); k += i; goto L_2020; } mess( mactop, (char*)mtxtac,8, &Ispec[k] ); smess( macttx, (char*)mtxtaa,210, idat, fdat ); Ispec[1] = -Macter[3]; return; } /* end of function */