/* ma27s.f -- translated by f2c (version 20200916). You must link the resulting object file with libf2c: on Microsoft Windows system, link with libf2c.lib; on Linux or Unix systems, link with .../path/to/libf2c.a -lm or, if you install libf2c.a in a standard place, with -lf2c -lm -- in that order, at the end of the command line, as in cc *.o -lf2c -lm Source for libf2c is in /netlib/f2c/libf2c.zip, e.g., http://www.netlib.org/f2c/libf2c.zip */ #include "f2c.h" /* Table of constant values */ static integer c__1 = 1; static integer c__2 = 2; /* COPYRIGHT (c) 1982 AEA Technology */ /* ######DATE 20 September 2001 */ /* September 2001: threadsafe version of MA27 */ /* 19/3/03. Array ICNTL in MA27G made assumed size. */ /* Subroutine */ int ma27i_(integer *icntl, real *cntl) { /* Stream number for error messages */ /* Parameter adjustments */ --cntl; --icntl; /* Function Body */ icntl[1] = 6; /* Stream number for diagnostic messages */ icntl[2] = 6; /* Control the level of diagnostic printing. */ /* 0 no printing */ /* 1 printing of scalar parameters and first parts of arrays. */ /* 2 printing of scalar parameters and whole of arrays. */ icntl[3] = 0; /* The largest integer such that all integers I in the range */ /* -ICNTL(4).LE.I.LE.ICNTL(4) can be handled by the shortest integer */ /* type in use. */ icntl[4] = 2139062143; /* Minimum number of eliminations in a step that is automatically */ /* accepted. if two adjacent steps can be combined and each has less */ /* eliminations then they are combined. */ icntl[5] = 1; /* Control whether direct or indirect access is used by MA27C/CD. */ /* Indirect access is employed in forward and back substitution */ /* respectively if the size of a block is less than */ /* ICNTL(IFRLVL+MIN(10,NPIV)) and ICNTL(IFRLVL+10+MIN(10,NPIV)) */ /* respectively, where NPIV is the number of pivots in the block. */ icntl[6] = 32639; icntl[7] = 32639; icntl[8] = 32639; icntl[9] = 32639; icntl[10] = 14; icntl[11] = 9; icntl[12] = 8; icntl[13] = 8; icntl[14] = 9; icntl[15] = 10; icntl[16] = 32639; icntl[17] = 32639; icntl[18] = 32639; icntl[19] = 32689; icntl[20] = 24; icntl[21] = 11; icntl[22] = 9; icntl[23] = 8; icntl[24] = 9; icntl[25] = 10; /* Not used */ icntl[26] = 0; icntl[27] = 0; icntl[28] = 0; icntl[29] = 0; icntl[30] = 0; /* Control threshold pivoting. */ cntl[1] = .1f; /* If a column of the reduced matrix has relative density greater than */ /* CNTL(2), it is forced into the root. All such columns are taken to */ /* have sparsity pattern equal to their merged patterns, so the fill */ /* and operation counts may be overestimated. */ cntl[2] = 1.f; /* An entry with absolute value less than CNTL(3) is never accepted as */ /* a 1x1 pivot or as the off-diagonal of a 2x2 pivot. */ cntl[3] = 0.f; /* Not used */ cntl[4] = 0.f; cntl[5] = 0.f; return 0; } /* ma27i_ */ /* Subroutine */ int ma27a_(integer *n, integer *nz, integer *irn, integer * icn, integer *iw, integer *liw, integer *ikeep, integer *iw1, integer *nsteps, integer *iflag, integer *icntl, real *cntl, integer *info, real *ops) { /* Format strings */ static char fmt_10[] = "(/,/,\002 ENTERING MA27A WITH N NZ " " LIW IFLAG\002,/,21x,i7,i7,i9,i7)"; static char fmt_20[] = "(\002 MATRIX NON-ZEROS\002,/,4(i9,i6),/,(i9,i6,i" "9,i6,i9,i6,i9,i6))"; static char fmt_30[] = "(\002 IKEEP(.,1)=\002,10i6,/,(12x,10i6))"; static char fmt_80[] = "(\002 **** ERROR RETURN FROM MA27A **** INFO(" "1)=\002,i3)"; static char fmt_90[] = "(\002 VALUE OF N OUT OF RANGE ... =\002,i10)"; static char fmt_110[] = "(\002 VALUE OF NZ OUT OF RANGE .. =\002,i10)"; static char fmt_150[] = "(\002 LIW TOO SMALL, MUST BE INCREASED FROM\002" ",i10,\002 TO AT LEAST\002,i10)"; static char fmt_170[] = "(/,\002 LEAVING MA27A WITH NSTEPS INFO(1) " "OPS IERROR\002,\002 NRLTOT NIRTOT\002,/,20x,2i7,f7.0,3i7,/,20x" ",\002 NRLNEC NIRNEC NRLADU NIRADU NCMPA\002,/,20x,6i7)"; static char fmt_180[] = "(\002 IKEEP(.,2)=\002,10i6,/,(12x,10i6))"; static char fmt_190[] = "(\002 IKEEP(.,3)=\002,10i6,/,(12x,10i6))"; /* System generated locals */ integer ikeep_dim1, ikeep_offset, iw1_dim1, iw1_offset, i__1; /* Builtin functions */ integer s_wsfe(cilist *), do_fio(integer *, char *, ftnlen), e_wsfe(void); /* Local variables */ static integer i__, k, l1, l2; extern /* Subroutine */ int ma27g_(integer *, integer *, integer *, integer *, integer *, integer *, integer *, integer *, integer *, integer *, integer *, integer *), ma27h_(integer *, integer *, integer *, integer *, integer *, integer *, integer *, integer *, integer *, integer *, integer *, integer *, real *), ma27j_( integer *, integer *, integer *, integer *, integer *, integer *, integer *, integer *, integer *, integer *, integer *, integer *, integer *), ma27k_(integer *, integer *, integer *, integer *, integer *, integer *, integer *, integer *, integer *, integer *), ma27l_(integer *, integer *, integer *, integer *, integer *, integer *, integer *, integer *, integer *), ma27m_(integer *, integer *, integer *, integer *, integer *, integer *, integer *, integer *, integer *, integer *, integer *, integer *, integer *, real *); static integer iwfr, lliw; /* Fortran I/O blocks */ static cilist io___2 = { 0, 0, 0, fmt_10, 0 }; static cilist io___4 = { 0, 0, 0, fmt_20, 0 }; static cilist io___5 = { 0, 0, 0, fmt_30, 0 }; static cilist io___10 = { 0, 0, 0, fmt_80, 0 }; static cilist io___11 = { 0, 0, 0, fmt_90, 0 }; static cilist io___12 = { 0, 0, 0, fmt_80, 0 }; static cilist io___13 = { 0, 0, 0, fmt_110, 0 }; static cilist io___14 = { 0, 0, 0, fmt_80, 0 }; static cilist io___15 = { 0, 0, 0, fmt_150, 0 }; static cilist io___16 = { 0, 0, 0, fmt_170, 0 }; static cilist io___17 = { 0, 0, 0, fmt_30, 0 }; static cilist io___18 = { 0, 0, 0, fmt_180, 0 }; static cilist io___19 = { 0, 0, 0, fmt_190, 0 }; /* THIS SUBROUTINE COMPUTES A MINIMUM DEGREE ORDERING OR ACCEPTS A GIVEN */ /* ORDERING. IT COMPUTES ASSOCIATED ASSEMBLY AND ELIMINATION */ /* INFORMATION FOR MA27B/BD. */ /* N MUST BE SET TO THE MATRIX ORDER. IT IS NOT ALTERED. */ /* NZ MUST BE SET TO THE NUMBER OF NON-ZEROS INPUT. IT IS NOT */ /* ALTERED. */ /* IRN(I),I=1,2,...,NZ MUST BE SET TO THE ROW NUMBERS OF THE */ /* NON-ZEROS. IT IS NOT ALTERED UNLESS IT IS EQUIVALENCED */ /* TO IW (SEE DESCRIPTION OF IW). */ /* ICN(I),I=1,2,...,NZ MUST BE SET TO THE COLUMN NUMBERS OF THE */ /* NON-ZEROS. IT IS NOT ALTERED UNLESS IT IS EQUIVALENCED */ /* TO IW (SEE DESCRIPTION OF IW). */ /* IW NEED NOT BE SET ON INPUT. IT IS USED FOR WORKSPACE. */ /* IRN(1) MAY BE EQUIVALENCED TO IW(1) AND ICN(1) MAY BE */ /* EQUIVALENCED TO IW(K), WHERE K.GT.NZ. */ /* LIW MUST BE SET TO THE LENGTH OF IW. IT MUST BE AT LEAST 2*NZ+3*N */ /* FOR THE IFLAG=0 ENTRY AND AT LEAST NZ+3*N FOR THE IFLAG=1 */ /* ENTRY. IT IS NOT ALTERED. */ /* IKEEP NEED NOT BE SET UNLESS AN ORDERING IS GIVEN, IN WHICH CASE */ /* IKEEP(I,1) MUST BE SET TO THE POSITION OF VARIABLE I IN THE */ /* ORDER. ON OUTPUT IKEEP CONTAINS INFORMATION NEEDED BY MA27B/BD. */ /* IKEEP(I,1) HOLDS THE POSITION OF VARIABLE I IN THE PIVOT ORDER. */ /* IKEEP(I,2), IKEEP(I,3) HOLD THE NUMBER OF ELIMINATIONS, ASSEMBLIES */ /* AT MAJOR STEP I, I=1,2,...,NSTEPS. NOTE THAT WHEN AN ORDER IS */ /* GIVEN IT MAY BE REPLACED BY ANOTHER ORDER THAT GIVES IDENTICAL */ /* NUMERICAL RESULTS. */ /* IW1 IS USED FOR WORKSPACE. */ /* NSTEPS NEED NOT BE SET. ON OUTPUT IT CONTAINS THE NUMBER OF MAJOR */ /* STEPS NEEDED FOR A DEFINITE MATRIX AND MUST BE PASSED UNCHANGED */ /* TO MA27B/BD. */ /* IFLAG MUST SET TO ZERO IF THE USER WANTS THE PIVOT ORDER CHOSEN */ /* AUTOMATICALLY AND TO ONE IF HE WANTS TO SPECIFY IT IN IKEEP. */ /* ICNTL is an INTEGER array of length 30 containing user options */ /* with integer values (defaults set in MA27I/ID) */ /* ICNTL(1) (LP) MUST BE SET TO THE STREAM NUMBER FOR ERROR MESSAGES. */ /* ERROR MESSAGES ARE SUPPRESSED IF ICNTL(1) IS NOT POSITIVE. */ /* IT IS NOT ALTERED. */ /* ICNTL(2) (MP) MUST BE SET TO THE STREAM NUMBER FOR DIAGNOSTIC */ /* MESSAGES. DIAGNOSTIC MESSAGES ARE SUPPRESSED IF ICNTL(2) IS NOT */ /* POSITIVE. IT IS NOT ALTERED. */ /* ICNTL(3) (LDIAG) CONTROLS THE LEVEL OF DIAGNOSTIC PRINTING. */ /* 0 NO PRINTING */ /* 1 PRINTING OF SCALAR PARAMETERS AND FIRST PARTS OF ARRAYS. */ /* 2 PRINTING OF SCALAR PARAMETERS AND WHOLE OF ARRAYS. */ /* ICNTL(4) (IOVFLO) IS THE LARGEST INTEGER SUCH THAT ALL INTEGERS */ /* I IN THE RANGE -IOVFLO.LE.I.LE.IOVFLO CAN BE HANDLED BY THE */ /* SHORTEST INTEGER TYPE IN USE. */ /* ICNT(5) (NEMIN) MUST BE SET TO THE MINIMUM NUMBER OF ELIMINATIONS */ /* IN A STEP THAT IS AUTOMATICALLY ACCEPTED. IF TWO ADJACENT STEPS */ /* CAN BE COMBINED AND EACH HAS LESS ELIMINATIONS THEN THEY ARE */ /* COMBINED. */ /* ICNTL(IFRLVL+I) I=1,20, (IFRLVL) MUST BE SET TO CONTROL WHETHER */ /* DIRECT OR INDIRECT ACCESS IS USED BY MA27C/CD. INDIRECT ACCESS */ /* IS EMPLOYED IN FORWARD AND BACK SUBSTITUTION RESPECTIVELY IF THE */ /* SIZE OF A BLOCK IS LESS THAN ICNTL(IFRLVL+(MIN(10,NPIV)) AND */ /* ICNTL(IFRLVL+10+MIN(10,NPIV)) RESPECTIVELY, WHERE NPIV IS THE */ /* NUMBER OF PIVOTS IN THE BLOCK. */ /* ICNTL(I) I=26,30 are not used. */ /* CNTL is an REAL array of length 5 containing user options */ /* with real values (defaults set in MA27I/ID) */ /* CNTL(1) (U) IS USED TO CONTROL THRESHOLD PIVOTING. IT IS NOT */ /* ALTERED. */ /* CNTL(2) (FRATIO) has default value 1.0. If a column of the */ /* reduced matrix has relative density greater than CNTL(2), it */ /* is forced into the root. All such columns are taken to have */ /* sparsity pattern equal to their merged patterns, so the fill */ /* and operation counts may be overestimated. */ /* CNTL(3) (PIVTOL) has default value 0.0. An entry with absolute */ /* value less than CNTL(3) is never accepted as a 1x1 pivot or */ /* as the off-diagonal of a 2x2 pivot. */ /* CNTL(I) I=4,5 are not used. */ /* INFO is an INTEGER array of length 20 which is used to return */ /* information to the user. */ /* INFO(1) (IFLAG) is an error return code, zero for success, greater */ /* than zero for a warning and less than zero for errors, see */ /* INFO(2). */ /* INFO(2) (IERROR) HOLDS ADDITIONAL INFORMATION IN THE EVENT OF ERRORS. */ /* IF INFO(1)=-3 INFO(2) HOLDS A LENGTH THAT MAY SUFFICE FOR IW. */ /* IF INFO(1)=-4 INFO(2) HOLDS A LENGTH THAT MAY SUFFICE FOR A. */ /* IF INFO(1)=-5 INFO(2) IS SET TO THE PIVOT STEP AT WHICH SINGULARITY */ /* WAS DETECTED. */ /* IF INFO(1)=-6 INFO(2) IS SET TO THE PIVOT STEP AT WHICH A CHANGE OF */ /* PIVOT SIGN WAS FOUND. */ /* IF INFO(1)= 1 INFO(2) HOLDS THE NUMBER OF FAULTY ENTRIES. */ /* IF INFO(1)= 2 INFO(2) IS SET TO THE NUMBER OF SIGNS CHANGES IN */ /* THE PIVOTS. */ /* IF INFO(1)=3 INFO(2) IS SET TO THE RANK OF THE MATRIX. */ /* INFO(3) and INFO(4) (NRLTOT and NIRTOT) REAL AND INTEGER STRORAGE */ /* RESPECTIVELY REQUIRED FOR THE FACTORIZATION IF NO COMPRESSES ARE */ /* ALLOWED. */ /* INFO(5) and INFO(6) (NRLNEC and NIRNEC) REAL AND INTEGER STORAGE */ /* RESPECTIVELY REQUIRED FOR THE FACTORIZATION IF COMPRESSES ARE */ /* ALLOWED AND THE MATRIX IS DEFINITE. */ /* INFO(7) and INFO(8) (NRLADU and NIRADU) REAL AND INTEGER STORAGE */ /* RESPECTIVELY FOR THE MATRIX FACTORS AS CALCULATED BY MA27A/AD */ /* FOR THE DEFINITE CASE. */ /* INFO(9) and INFO(10) (NRLBDU and NIRBDU) REAL AND INTEGER STORAGE */ /* RESPECTIVELY FOR THE MATRIX FACTORS AS FOUND BY MA27B/BD. */ /* INFO(11) (NCMPA) ACCUMULATES THE NUMBER OF TIMES THE ARRAY IW IS */ /* COMPRESSED BY MA27A/AD. */ /* INFO(12) and INFO(13) (NCMPBR and NCMPBI) ACCUMULATE THE NUMBER */ /* OF COMPRESSES OF THE REALS AND INTEGERS PERFORMED BY MA27B/BD. */ /* INFO(14) (NTWO) IS USED BY MA27B/BD TO RECORD THE NUMBER OF 2*2 */ /* PIVOTS USED. */ /* INFO(15) (NEIG) IS USED BY ME27B/BD TO RECORD THE NUMBER OF */ /* NEGATIVE EIGENVALUES OF A. */ /* INFO(16) to INFO(20) are not used. */ /* OPS ACCUMULATES THE NO. OF MULTIPLY/ADD PAIRS NEEDED TO CREATE THE */ /* TRIANGULAR FACTORIZATION, IN THE DEFINITE CASE. */ /* .. Scalar Arguments .. */ /* .. */ /* .. Array Arguments .. */ /* .. */ /* .. Local Scalars .. */ /* .. */ /* .. External Subroutines .. */ /* .. */ /* .. Intrinsic Functions .. */ /* .. */ /* .. Executable Statements .. */ /* Parameter adjustments */ iw1_dim1 = *n; iw1_offset = 1 + iw1_dim1; iw1 -= iw1_offset; ikeep_dim1 = *n; ikeep_offset = 1 + ikeep_dim1; ikeep -= ikeep_offset; --irn; --icn; --iw; --icntl; --cntl; --info; /* Function Body */ for (i__ = 1; i__ <= 15; ++i__) { info[i__] = 0; /* L5: */ } if (icntl[3] <= 0 || icntl[2] <= 0) { goto L40; } /* PRINT INPUT VARIABLES. */ io___2.ciunit = icntl[2]; s_wsfe(&io___2); do_fio(&c__1, (char *)&(*n), (ftnlen)sizeof(integer)); do_fio(&c__1, (char *)&(*nz), (ftnlen)sizeof(integer)); do_fio(&c__1, (char *)&(*liw), (ftnlen)sizeof(integer)); do_fio(&c__1, (char *)&(*iflag), (ftnlen)sizeof(integer)); e_wsfe(); *nsteps = 0; k = min(8,*nz); if (icntl[3] > 1) { k = *nz; } if (k > 0) { io___4.ciunit = icntl[2]; s_wsfe(&io___4); i__1 = k; for (i__ = 1; i__ <= i__1; ++i__) { do_fio(&c__1, (char *)&irn[i__], (ftnlen)sizeof(integer)); do_fio(&c__1, (char *)&icn[i__], (ftnlen)sizeof(integer)); } e_wsfe(); } k = min(10,*n); if (icntl[3] > 1) { k = *n; } if (*iflag == 1 && k > 0) { io___5.ciunit = icntl[2]; s_wsfe(&io___5); i__1 = k; for (i__ = 1; i__ <= i__1; ++i__) { do_fio(&c__1, (char *)&ikeep[i__ + ikeep_dim1], (ftnlen)sizeof( integer)); } e_wsfe(); } L40: if (*n < 1 || *n > icntl[4]) { goto L70; } if (*nz < 0) { goto L100; } lliw = *liw - (*n << 1); l1 = lliw + 1; l2 = l1 + *n; if (*iflag == 1) { goto L50; } if (*liw < (*nz << 1) + *n * 3 + 1) { goto L130; } /* SORT */ ma27g_(n, nz, &irn[1], &icn[1], &iw[1], &lliw, &iw1[iw1_offset], &iw1[( iw1_dim1 << 1) + 1], &iw[l1], &iwfr, &icntl[1], &info[1]); /* ANALYZE USING MINIMUM DEGREE ORDERING */ ma27h_(n, &iw1[iw1_offset], &iw[1], &lliw, &iwfr, &iw[l1], &iw[l2], & ikeep[(ikeep_dim1 << 1) + 1], &ikeep[ikeep_dim1 * 3 + 1], &ikeep[ ikeep_offset], &icntl[4], &info[11], &cntl[2]); goto L60; /* SORT USING GIVEN ORDER */ L50: if (*liw < *nz + *n * 3 + 1) { goto L120; } ma27j_(n, nz, &irn[1], &icn[1], &ikeep[ikeep_offset], &iw[1], &lliw, &iw1[ iw1_offset], &iw1[(iw1_dim1 << 1) + 1], &iw[l1], &iwfr, &icntl[1], &info[1]); /* ANALYZE USING GIVEN ORDER */ ma27k_(n, &iw1[iw1_offset], &iw[1], &lliw, &iwfr, &ikeep[ikeep_offset], & ikeep[(ikeep_dim1 << 1) + 1], &iw[l1], &iw[l2], &info[11]); /* PERFORM DEPTH-FIRST SEARCH OF ASSEMBLY TREE */ L60: ma27l_(n, &iw1[iw1_offset], &iw[l1], &ikeep[ikeep_offset], &ikeep[( ikeep_dim1 << 1) + 1], &ikeep[ikeep_dim1 * 3 + 1], &iw[l2], nsteps, &icntl[5]); /* EVALUATE STORAGE AND OPERATION COUNT REQUIRED BY MA27B/BD IN THE */ /* DEFINITE CASE. */ /* SET IW(1) SO THAT ARRAYS IW AND IRN CAN BE TESTED FOR EQUIVALENCE */ /* IN MA27M/MD. */ if (*nz >= 1) { iw[1] = irn[1] + 1; } ma27m_(n, nz, &irn[1], &icn[1], &ikeep[ikeep_offset], &ikeep[ikeep_dim1 * 3 + 1], &ikeep[(ikeep_dim1 << 1) + 1], &iw[l2], nsteps, &iw1[ iw1_offset], &iw1[(iw1_dim1 << 1) + 1], &iw[1], &info[1], ops); goto L160; L70: info[1] = -1; if (icntl[1] > 0) { io___10.ciunit = icntl[1]; s_wsfe(&io___10); do_fio(&c__1, (char *)&info[1], (ftnlen)sizeof(integer)); e_wsfe(); } if (icntl[1] > 0) { io___11.ciunit = icntl[1]; s_wsfe(&io___11); do_fio(&c__1, (char *)&(*n), (ftnlen)sizeof(integer)); e_wsfe(); } goto L160; L100: info[1] = -2; if (icntl[1] > 0) { io___12.ciunit = icntl[1]; s_wsfe(&io___12); do_fio(&c__1, (char *)&info[1], (ftnlen)sizeof(integer)); e_wsfe(); } if (icntl[1] > 0) { io___13.ciunit = icntl[1]; s_wsfe(&io___13); do_fio(&c__1, (char *)&(*nz), (ftnlen)sizeof(integer)); e_wsfe(); } goto L160; L120: info[2] = *nz + *n * 3 + 1; goto L140; L130: info[2] = (*nz << 1) + *n * 3 + 1; L140: info[1] = -3; if (icntl[1] > 0) { io___14.ciunit = icntl[1]; s_wsfe(&io___14); do_fio(&c__1, (char *)&info[1], (ftnlen)sizeof(integer)); e_wsfe(); } if (icntl[1] > 0) { io___15.ciunit = icntl[1]; s_wsfe(&io___15); do_fio(&c__1, (char *)&(*liw), (ftnlen)sizeof(integer)); do_fio(&c__1, (char *)&info[2], (ftnlen)sizeof(integer)); e_wsfe(); } L160: if (icntl[3] <= 0 || icntl[2] <= 0 || info[1] < 0) { goto L200; } /* PRINT PARAMETER VALUES ON EXIT. */ io___16.ciunit = icntl[2]; s_wsfe(&io___16); do_fio(&c__1, (char *)&(*nsteps), (ftnlen)sizeof(integer)); do_fio(&c__1, (char *)&info[1], (ftnlen)sizeof(integer)); do_fio(&c__1, (char *)&(*ops), (ftnlen)sizeof(real)); do_fio(&c__1, (char *)&info[2], (ftnlen)sizeof(integer)); do_fio(&c__1, (char *)&info[3], (ftnlen)sizeof(integer)); do_fio(&c__1, (char *)&info[4], (ftnlen)sizeof(integer)); do_fio(&c__1, (char *)&info[5], (ftnlen)sizeof(integer)); do_fio(&c__1, (char *)&info[6], (ftnlen)sizeof(integer)); do_fio(&c__1, (char *)&info[7], (ftnlen)sizeof(integer)); do_fio(&c__1, (char *)&info[8], (ftnlen)sizeof(integer)); do_fio(&c__1, (char *)&info[11], (ftnlen)sizeof(integer)); e_wsfe(); k = min(9,*n); if (icntl[3] > 1) { k = *n; } if (k > 0) { io___17.ciunit = icntl[2]; s_wsfe(&io___17); i__1 = k; for (i__ = 1; i__ <= i__1; ++i__) { do_fio(&c__1, (char *)&ikeep[i__ + ikeep_dim1], (ftnlen)sizeof( integer)); } e_wsfe(); } k = min(k,*nsteps); if (k > 0) { io___18.ciunit = icntl[2]; s_wsfe(&io___18); i__1 = k; for (i__ = 1; i__ <= i__1; ++i__) { do_fio(&c__1, (char *)&ikeep[i__ + (ikeep_dim1 << 1)], (ftnlen) sizeof(integer)); } e_wsfe(); } if (k > 0) { io___19.ciunit = icntl[2]; s_wsfe(&io___19); i__1 = k; for (i__ = 1; i__ <= i__1; ++i__) { do_fio(&c__1, (char *)&ikeep[i__ + ikeep_dim1 * 3], (ftnlen) sizeof(integer)); } e_wsfe(); } L200: return 0; } /* ma27a_ */ /* Subroutine */ int ma27b_(integer *n, integer *nz, integer *irn, integer * icn, real *a, integer *la, integer *iw, integer *liw, integer *ikeep, integer *nsteps, integer *maxfrt, integer *iw1, integer *icntl, real * cntl, integer *info) { /* Format strings */ static char fmt_10[] = "(/,/,\002 ENTERING MA27B WITH N NZ " " LA LIW\002,\002 NSTEPS U\002,/,21x,i7,i7,i9,i9,i7,1" "pe10.2)"; static char fmt_20[] = "(\002 MATRIX NON-ZEROS\002,/,1x,2(1p,e16.3,2i6)," "/,(1x,1p,e16.3,2i6,1p,e16.3,2i6))"; static char fmt_30[] = "(\002 IKEEP(.,1)=\002,10i6,/,(12x,10i6))"; static char fmt_40[] = "(\002 IKEEP(.,2)=\002,10i6,/,(12x,10i6))"; static char fmt_50[] = "(\002 IKEEP(.,3)=\002,10i6,/,(12x,10i6))"; static char fmt_65[] = "(\002 *** WARNING MESSAGE FROM SUBROUTINE MA27" "B\002,\002 *** INFO(1) =\002,i2,/,5x,\002MATRIX IS SINGULAR. RA" "NK=\002,i5)"; static char fmt_80[] = "(\002 **** ERROR RETURN FROM MA27B **** INFO(" "1)=\002,i3)"; static char fmt_90[] = "(\002 VALUE OF N OUT OF RANGE ... =\002,i10)"; static char fmt_110[] = "(\002 VALUE OF NZ OUT OF RANGE .. =\002,i10)"; static char fmt_140[] = "(\002 LIW TOO SMALL, MUST BE INCREASED FROM\002" ",i10,\002 TO\002,\002 AT LEAST\002,i10)"; static char fmt_170[] = "(\002 LA TOO SMALL, MUST BE INCREASED FROM \002" ",i10,\002 TO\002,\002 AT LEAST\002,i10)"; static char fmt_190[] = "(\002 ZERO PIVOT AT STAGE\002,i10,\002 WHEN INP" "UT MATRIX DECLARED DEFINITE\002)"; static char fmt_210[] = "(\002 CHANGE IN SIGN OF PIVOT ENCOUNTERED\002" ",\002 WHEN FACTORING ALLEGEDLY DEFINITE MATRIX\002)"; static char fmt_230[] = "(/,\002 LEAVING MA27B WITH\002,/,10x,\002 MAX" "FRT INFO(1) NRLBDU NIRBDU NCMPBR\002,\002 NCMPBI NTWO IERRO" "R\002,/,11x,8i7)"; static char fmt_250[] = "(\002 BLOCK PIVOT =\002,i8,\002 NROWS =\002,i8" ",\002 NCOLS =\002,i8)"; static char fmt_260[] = "(\002 COLUMN INDICES =\002,10i6,/,(17x,10i6))"; static char fmt_270[] = "(\002 REAL ENTRIES .. EACH ROW STARTS ON A NEW " "LINE\002)"; static char fmt_280[] = "(1p,5e13.3)"; /* System generated locals */ integer ikeep_dim1, ikeep_offset, i__1, i__2, i__3; /* Builtin functions */ integer s_wsfe(cilist *), do_fio(integer *, char *, ftnlen), e_wsfe(void); /* Local variables */ static integer i__, k, j1, j2, jj, kz, nz1, len, iblk, kblk; extern /* Subroutine */ int ma27n_(integer *, integer *, integer *, real * , integer *, integer *, integer *, integer *, integer *, integer * , integer *, integer *, integer *), ma27o_(integer *, integer *, real *, integer *, integer *, integer *, integer *, integer *, integer *, integer *, integer *, integer *, integer *, real *, integer *); static integer ipos, iapos, ncols, irows, nrows; /* Fortran I/O blocks */ static cilist io___20 = { 0, 0, 0, fmt_10, 0 }; static cilist io___22 = { 0, 0, 0, fmt_20, 0 }; static cilist io___24 = { 0, 0, 0, fmt_30, 0 }; static cilist io___26 = { 0, 0, 0, fmt_40, 0 }; static cilist io___27 = { 0, 0, 0, fmt_50, 0 }; static cilist io___29 = { 0, 0, 0, fmt_65, 0 }; static cilist io___30 = { 0, 0, 0, fmt_80, 0 }; static cilist io___31 = { 0, 0, 0, fmt_90, 0 }; static cilist io___32 = { 0, 0, 0, fmt_80, 0 }; static cilist io___33 = { 0, 0, 0, fmt_110, 0 }; static cilist io___34 = { 0, 0, 0, fmt_80, 0 }; static cilist io___35 = { 0, 0, 0, fmt_140, 0 }; static cilist io___36 = { 0, 0, 0, fmt_80, 0 }; static cilist io___37 = { 0, 0, 0, fmt_170, 0 }; static cilist io___38 = { 0, 0, 0, fmt_80, 0 }; static cilist io___39 = { 0, 0, 0, "(A)", 0 }; static cilist io___40 = { 0, 0, 0, fmt_80, 0 }; static cilist io___41 = { 0, 0, 0, fmt_190, 0 }; static cilist io___42 = { 0, 0, 0, fmt_80, 0 }; static cilist io___43 = { 0, 0, 0, fmt_210, 0 }; static cilist io___44 = { 0, 0, 0, fmt_230, 0 }; static cilist io___52 = { 0, 0, 0, fmt_250, 0 }; static cilist io___54 = { 0, 0, 0, fmt_260, 0 }; static cilist io___56 = { 0, 0, 0, fmt_270, 0 }; static cilist io___59 = { 0, 0, 0, fmt_280, 0 }; /* THIS SUBROUTINE COMPUTES THE FACTORISATION OF THE MATRIX INPUT IN */ /* A,IRN,ICN USING INFORMATION (IN IKEEP) FROM MA27A/AD. */ /* N MUST BE SET TO THE ORDER OF THE MATRIX. IT IS NOT ALTERED. */ /* NZ MUST BE SET TO THE NUMBER OF NON-ZEROS INPUT. IT IS NOT */ /* ALTERED. */ /* IRN,ICN,A. ENTRY K (K=1,NZ) OF IRN,ICN MUST BE SET TO THE ROW */ /* AND COLUMN INDEX RESPECTIVELY OF THE NON-ZERO IN A. */ /* IRN AND ICN ARE UNALTERED BY MA27B/BD. */ /* ON EXIT, ENTRIES 1 TO NRLBDU OF A HOLD REAL INFORMATION */ /* ON THE FACTORS AND SHOULD BE PASSED UNCHANGED TO MA27C/CD. */ /* LA LENGTH OF ARRAY A. AN INDICATION OF A SUITABLE VALUE, */ /* SUFFICIENT FOR FACTORIZATION OF A DEFINITE MATRIX, WILL */ /* HAVE BEEN PROVIDED IN NRLNEC AND NRLTOT BY MA27A/AD. */ /* IT IS NOT ALTERED BY MA27B/BD. */ /* IW NEED NOT BE SET ON ENTRY. USED AS A WORK ARRAY BY MA27B/BD. */ /* ON EXIT, ENTRIES 1 TO NIRBDU HOLD INTEGER INFORMATION ON THE */ /* FACTORS AND SHOULD BE PASSED UNCHANGED TO MA27C/CD. */ /* LIW LENGTH OF ARRAY IW. AN INDICATION OF A SUITABLE VALUE WILL */ /* HAVE BEEN PROVIDED IN NIRNEC AND NIRTOT BY MA27A/AD. */ /* IT IS NOT ALTERED BY MA27B/BD. */ /* IKEEP MUST BE UNCHANGED SINCE THE CALL TO MA27A/AD. */ /* IT IS NOT ALTERED BY MA27B/BD. */ /* NSTEPS MUST BE UNCHANGED SINCE THE CALL TO MA27A/AD. */ /* IT IS NOT ALTERED BY MA27B/BD. */ /* MAXFRT NEED NOT BE SET AND MUST BE PASSED UNCHANGED TO MA27C/CD. */ /* IT IS THE MAXIMUM SIZE OF THE FRONTAL MATRIX GENERATED DURING */ /* THE DECOMPOSITION. */ /* IW1 USED AS WORKSPACE BY MA27B/BD. */ /* ICNTL is an INTEGER array of length 30, see MA27A/AD. */ /* CNTL is a REAL array of length 5, see MA27A/AD. */ /* INFO is an INTEGER array of length 20, see MA27A/AD. */ /* .. Scalar Arguments .. */ /* .. */ /* .. Array Arguments .. */ /* .. */ /* .. Local Scalars .. */ /* .. */ /* .. External Subroutines .. */ /* .. */ /* .. Intrinsic Functions .. */ /* .. */ /* .. Executable Statements .. */ /* Parameter adjustments */ --iw1; ikeep_dim1 = *n; ikeep_offset = 1 + ikeep_dim1; ikeep -= ikeep_offset; --irn; --icn; --a; --iw; --icntl; --cntl; --info; /* Function Body */ info[1] = 0; if (icntl[3] <= 0 || icntl[2] <= 0) { goto L60; } /* PRINT INPUT PARAMETERS. */ io___20.ciunit = icntl[2]; s_wsfe(&io___20); do_fio(&c__1, (char *)&(*n), (ftnlen)sizeof(integer)); do_fio(&c__1, (char *)&(*nz), (ftnlen)sizeof(integer)); do_fio(&c__1, (char *)&(*la), (ftnlen)sizeof(integer)); do_fio(&c__1, (char *)&(*liw), (ftnlen)sizeof(integer)); do_fio(&c__1, (char *)&(*nsteps), (ftnlen)sizeof(integer)); do_fio(&c__1, (char *)&cntl[1], (ftnlen)sizeof(real)); e_wsfe(); kz = min(6,*nz); if (icntl[3] > 1) { kz = *nz; } if (*nz > 0) { io___22.ciunit = icntl[2]; s_wsfe(&io___22); i__1 = kz; for (k = 1; k <= i__1; ++k) { do_fio(&c__1, (char *)&a[k], (ftnlen)sizeof(real)); do_fio(&c__1, (char *)&irn[k], (ftnlen)sizeof(integer)); do_fio(&c__1, (char *)&icn[k], (ftnlen)sizeof(integer)); } e_wsfe(); } k = min(9,*n); if (icntl[3] > 1) { k = *n; } if (k > 0) { io___24.ciunit = icntl[2]; s_wsfe(&io___24); i__1 = k; for (i__ = 1; i__ <= i__1; ++i__) { do_fio(&c__1, (char *)&ikeep[i__ + ikeep_dim1], (ftnlen)sizeof( integer)); } e_wsfe(); } k = min(k,*nsteps); if (k > 0) { io___26.ciunit = icntl[2]; s_wsfe(&io___26); i__1 = k; for (i__ = 1; i__ <= i__1; ++i__) { do_fio(&c__1, (char *)&ikeep[i__ + (ikeep_dim1 << 1)], (ftnlen) sizeof(integer)); } e_wsfe(); } if (k > 0) { io___27.ciunit = icntl[2]; s_wsfe(&io___27); i__1 = k; for (i__ = 1; i__ <= i__1; ++i__) { do_fio(&c__1, (char *)&ikeep[i__ + ikeep_dim1 * 3], (ftnlen) sizeof(integer)); } e_wsfe(); } L60: if (*n < 1 || *n > icntl[4]) { goto L70; } if (*nz < 0) { goto L100; } if (*liw < *nz) { goto L120; } if (*la < *nz + *n) { goto L150; } if (*nsteps < 1 || *nsteps > *n) { goto L175; } /* SORT */ ma27n_(n, nz, &nz1, &a[1], la, &irn[1], &icn[1], &iw[1], liw, &ikeep[ ikeep_offset], &iw1[1], &icntl[1], &info[1]); if (info[1] == -3) { goto L130; } if (info[1] == -4) { goto L160; } /* FACTORIZE */ ma27o_(n, &nz1, &a[1], la, &iw[1], liw, &ikeep[ikeep_offset], &ikeep[ ikeep_dim1 * 3 + 1], nsteps, maxfrt, &ikeep[(ikeep_dim1 << 1) + 1] , &iw1[1], &icntl[1], &cntl[1], &info[1]); if (info[1] == -3) { goto L130; } if (info[1] == -4) { goto L160; } if (info[1] == -5) { goto L180; } if (info[1] == -6) { goto L200; } /* **** WARNING MESSAGE **** */ if (info[1] == 3 && icntl[2] > 0) { io___29.ciunit = icntl[2]; s_wsfe(&io___29); do_fio(&c__1, (char *)&info[1], (ftnlen)sizeof(integer)); do_fio(&c__1, (char *)&info[2], (ftnlen)sizeof(integer)); e_wsfe(); } goto L220; /* **** ERROR RETURNS **** */ L70: info[1] = -1; if (icntl[1] > 0) { io___30.ciunit = icntl[1]; s_wsfe(&io___30); do_fio(&c__1, (char *)&info[1], (ftnlen)sizeof(integer)); e_wsfe(); } if (icntl[1] > 0) { io___31.ciunit = icntl[1]; s_wsfe(&io___31); do_fio(&c__1, (char *)&(*n), (ftnlen)sizeof(integer)); e_wsfe(); } goto L220; L100: info[1] = -2; if (icntl[1] > 0) { io___32.ciunit = icntl[1]; s_wsfe(&io___32); do_fio(&c__1, (char *)&info[1], (ftnlen)sizeof(integer)); e_wsfe(); } if (icntl[1] > 0) { io___33.ciunit = icntl[1]; s_wsfe(&io___33); do_fio(&c__1, (char *)&(*nz), (ftnlen)sizeof(integer)); e_wsfe(); } goto L220; L120: info[1] = -3; info[2] = *nz; L130: if (icntl[1] > 0) { io___34.ciunit = icntl[1]; s_wsfe(&io___34); do_fio(&c__1, (char *)&info[1], (ftnlen)sizeof(integer)); e_wsfe(); } if (icntl[1] > 0) { io___35.ciunit = icntl[1]; s_wsfe(&io___35); do_fio(&c__1, (char *)&(*liw), (ftnlen)sizeof(integer)); do_fio(&c__1, (char *)&info[2], (ftnlen)sizeof(integer)); e_wsfe(); } goto L220; L150: info[1] = -4; info[2] = *nz + *n; L160: if (icntl[1] > 0) { io___36.ciunit = icntl[1]; s_wsfe(&io___36); do_fio(&c__1, (char *)&info[1], (ftnlen)sizeof(integer)); e_wsfe(); } if (icntl[1] > 0) { io___37.ciunit = icntl[1]; s_wsfe(&io___37); do_fio(&c__1, (char *)&(*la), (ftnlen)sizeof(integer)); do_fio(&c__1, (char *)&info[2], (ftnlen)sizeof(integer)); e_wsfe(); } goto L220; L175: info[1] = -7; if (icntl[1] > 0) { io___38.ciunit = icntl[1]; s_wsfe(&io___38); do_fio(&c__1, (char *)&info[1], (ftnlen)sizeof(integer)); e_wsfe(); } if (icntl[1] > 0) { io___39.ciunit = icntl[1]; s_wsfe(&io___39); do_fio(&c__1, " NSTEPS is out of range", (ftnlen)23); e_wsfe(); } goto L220; L180: if (icntl[1] > 0) { io___40.ciunit = icntl[1]; s_wsfe(&io___40); do_fio(&c__1, (char *)&info[1], (ftnlen)sizeof(integer)); e_wsfe(); } if (icntl[1] > 0) { io___41.ciunit = icntl[1]; s_wsfe(&io___41); do_fio(&c__1, (char *)&info[2], (ftnlen)sizeof(integer)); e_wsfe(); } goto L220; L200: if (icntl[1] > 0) { io___42.ciunit = icntl[1]; s_wsfe(&io___42); do_fio(&c__1, (char *)&info[1], (ftnlen)sizeof(integer)); e_wsfe(); } if (icntl[1] > 0) { io___43.ciunit = icntl[1]; s_wsfe(&io___43); e_wsfe(); } L220: if (icntl[3] <= 0 || icntl[2] <= 0 || info[1] < 0) { goto L310; } /* PRINT OUTPUT PARAMETERS. */ io___44.ciunit = icntl[2]; s_wsfe(&io___44); do_fio(&c__1, (char *)&(*maxfrt), (ftnlen)sizeof(integer)); do_fio(&c__1, (char *)&info[1], (ftnlen)sizeof(integer)); do_fio(&c__1, (char *)&info[9], (ftnlen)sizeof(integer)); do_fio(&c__1, (char *)&info[10], (ftnlen)sizeof(integer)); do_fio(&c__1, (char *)&info[12], (ftnlen)sizeof(integer)); do_fio(&c__1, (char *)&info[13], (ftnlen)sizeof(integer)); do_fio(&c__1, (char *)&info[14], (ftnlen)sizeof(integer)); do_fio(&c__1, (char *)&info[2], (ftnlen)sizeof(integer)); e_wsfe(); if (info[1] < 0) { goto L310; } /* PRINT OUT MATRIX FACTORS FROM MA27B/BD. */ kblk = abs(iw[1]); if (kblk == 0) { goto L310; } if (icntl[3] == 1) { kblk = 1; } ipos = 2; iapos = 1; i__1 = kblk; for (iblk = 1; iblk <= i__1; ++iblk) { ncols = iw[ipos]; nrows = iw[ipos + 1]; j1 = ipos + 2; if (ncols > 0) { goto L240; } ncols = -ncols; nrows = 1; --j1; L240: io___52.ciunit = icntl[2]; s_wsfe(&io___52); do_fio(&c__1, (char *)&iblk, (ftnlen)sizeof(integer)); do_fio(&c__1, (char *)&nrows, (ftnlen)sizeof(integer)); do_fio(&c__1, (char *)&ncols, (ftnlen)sizeof(integer)); e_wsfe(); j2 = j1 + ncols - 1; ipos = j2 + 1; io___54.ciunit = icntl[2]; s_wsfe(&io___54); i__2 = j2; for (jj = j1; jj <= i__2; ++jj) { do_fio(&c__1, (char *)&iw[jj], (ftnlen)sizeof(integer)); } e_wsfe(); io___56.ciunit = icntl[2]; s_wsfe(&io___56); e_wsfe(); len = ncols; i__2 = nrows; for (irows = 1; irows <= i__2; ++irows) { j1 = iapos; j2 = iapos + len - 1; io___59.ciunit = icntl[2]; s_wsfe(&io___59); i__3 = j2; for (jj = j1; jj <= i__3; ++jj) { do_fio(&c__1, (char *)&a[jj], (ftnlen)sizeof(real)); } e_wsfe(); --len; iapos = j2 + 1; /* L290: */ } /* L300: */ } L310: return 0; } /* ma27b_ */ /* Subroutine */ int ma27c_(integer *n, real *a, integer *la, integer *iw, integer *liw, real *w, integer *maxfrt, real *rhs, integer *iw1, integer *nsteps, integer *icntl, integer *info) { /* Format strings */ static char fmt_10[] = "(/,/,\002 ENTERING MA27C WITH N LA " "LIW MAXFRT\002,\002 NSTEPS\002,/,21x,5i7)"; static char fmt_30[] = "(\002 BLOCK PIVOT =\002,i8,\002 NROWS =\002,i8" ",\002 NCOLS =\002,i8)"; static char fmt_40[] = "(\002 COLUMN INDICES =\002,10i6,/,(17x,10i6))"; static char fmt_50[] = "(\002 REAL ENTRIES .. EACH ROW STARTS ON A NEW L" "INE\002)"; static char fmt_60[] = "(1p,5e13.3)"; static char fmt_100[] = "(\002 RHS\002,1p,5e13.3,/,(4x,1p,5e13.3))"; static char fmt_160[] = "(/,/,\002 LEAVING MA27C WITH\002)"; /* System generated locals */ integer i__1, i__2, i__3; /* Builtin functions */ integer s_wsfe(cilist *), do_fio(integer *, char *, ftnlen), e_wsfe(void); /* Local variables */ static integer i__, k, j1, j2, jj, len, iblk, kblk, nblk; extern /* Subroutine */ int ma27q_(integer *, real *, integer *, integer * , integer *, real *, integer *, real *, integer *, integer *, integer *, integer *), ma27r_(integer *, real *, integer *, integer *, integer *, real *, integer *, real *, integer *, integer *, integer *, integer *); static integer ipos, iapos, ncols, latop, irows, nrows; /* Fortran I/O blocks */ static cilist io___60 = { 0, 0, 0, fmt_10, 0 }; static cilist io___68 = { 0, 0, 0, fmt_30, 0 }; static cilist io___70 = { 0, 0, 0, fmt_40, 0 }; static cilist io___72 = { 0, 0, 0, fmt_50, 0 }; static cilist io___75 = { 0, 0, 0, fmt_60, 0 }; static cilist io___77 = { 0, 0, 0, fmt_100, 0 }; static cilist io___81 = { 0, 0, 0, fmt_160, 0 }; static cilist io___82 = { 0, 0, 0, fmt_100, 0 }; /* THIS SUBROUTINE USES THE FACTORISATION OF THE MATRIX IN A,IW TO */ /* SOLVE A SYSTEM OF EQUATIONS. */ /* N MUST BE SET TO THE ORDER OF THE MATRIX. IT IS NOT ALTERED. */ /* A,IW HOLD INFORMATION ON THE FACTORS AND MUST BE UNCHANGED SINCE */ /* THE CALL TO MA27B/BD. THEY ARE NOT ALTERED BY MA27C/CDD. */ /* LA,LIW MUST BE SET TO THE LENGTHS OF A,IW RESPECTIVELY. THEY */ /* ARE NOT ALTERED. */ /* W USED AS A WORK ARRAY. */ /* MAXFRT IS THE LENGTH OF W AND MUST BE PASSED UNCHANGED FROM THE */ /* CALL TO MA27B/BD. IT IS NOT ALTERED. */ /* RHS MUST BE SET TO THE RIGHT HAND SIDE FOR THE EQUATIONS BEING */ /* SOLVED. ON EXIT, THIS ARRAY WILL HOLD THE SOLUTION. */ /* IW1 USED AS A WORK ARRAY. */ /* NSTEPS IS THE LENGTH OF IW1 WHICH MUST BE AT LEAST THE ABSOLUTE */ /* VALUE OF IW(1). IT IS NOT ALTERED. */ /* ICNTL is an INTEGER array of length 30, see MA27A/AD. */ /* INFO is an INTEGER array of length 20, see MA27A/AD. */ /* .. Scalar Arguments .. */ /* .. */ /* .. Array Arguments .. */ /* .. */ /* .. Local Scalars .. */ /* .. */ /* .. External Subroutines .. */ /* .. */ /* .. Intrinsic Functions .. */ /* .. */ /* .. Executable Statements .. */ /* Parameter adjustments */ --rhs; --a; --iw; --w; --iw1; --icntl; --info; /* Function Body */ info[1] = 0; if (icntl[3] <= 0 || icntl[2] <= 0) { goto L110; } /* PRINT INPUT PARAMETERS. */ io___60.ciunit = icntl[2]; s_wsfe(&io___60); do_fio(&c__1, (char *)&(*n), (ftnlen)sizeof(integer)); do_fio(&c__1, (char *)&(*la), (ftnlen)sizeof(integer)); do_fio(&c__1, (char *)&(*liw), (ftnlen)sizeof(integer)); do_fio(&c__1, (char *)&(*maxfrt), (ftnlen)sizeof(integer)); do_fio(&c__1, (char *)&(*nsteps), (ftnlen)sizeof(integer)); e_wsfe(); /* PRINT OUT MATRIX FACTORS FROM MA27B/BD. */ kblk = abs(iw[1]); if (kblk == 0) { goto L90; } if (icntl[3] == 1) { kblk = 1; } ipos = 2; iapos = 1; i__1 = kblk; for (iblk = 1; iblk <= i__1; ++iblk) { ncols = iw[ipos]; nrows = iw[ipos + 1]; j1 = ipos + 2; if (ncols > 0) { goto L20; } ncols = -ncols; nrows = 1; --j1; L20: io___68.ciunit = icntl[2]; s_wsfe(&io___68); do_fio(&c__1, (char *)&iblk, (ftnlen)sizeof(integer)); do_fio(&c__1, (char *)&nrows, (ftnlen)sizeof(integer)); do_fio(&c__1, (char *)&ncols, (ftnlen)sizeof(integer)); e_wsfe(); j2 = j1 + ncols - 1; ipos = j2 + 1; io___70.ciunit = icntl[2]; s_wsfe(&io___70); i__2 = j2; for (jj = j1; jj <= i__2; ++jj) { do_fio(&c__1, (char *)&iw[jj], (ftnlen)sizeof(integer)); } e_wsfe(); io___72.ciunit = icntl[2]; s_wsfe(&io___72); e_wsfe(); len = ncols; i__2 = nrows; for (irows = 1; irows <= i__2; ++irows) { j1 = iapos; j2 = iapos + len - 1; io___75.ciunit = icntl[2]; s_wsfe(&io___75); i__3 = j2; for (jj = j1; jj <= i__3; ++jj) { do_fio(&c__1, (char *)&a[jj], (ftnlen)sizeof(real)); } e_wsfe(); --len; iapos = j2 + 1; /* L70: */ } /* L80: */ } L90: k = min(10,*n); if (icntl[3] > 1) { k = *n; } if (*n > 0) { io___77.ciunit = icntl[2]; s_wsfe(&io___77); i__1 = k; for (i__ = 1; i__ <= i__1; ++i__) { do_fio(&c__1, (char *)&rhs[i__], (ftnlen)sizeof(real)); } e_wsfe(); } L110: if (iw[1] < 0) { goto L130; } nblk = iw[1]; if (nblk > 0) { goto L140; } /* WE HAVE A ZERO MATRIX */ i__1 = *n; for (i__ = 1; i__ <= i__1; ++i__) { rhs[i__] = 0.f; /* L120: */ } goto L150; L130: nblk = -iw[1]; /* FORWARD SUBSTITUTION */ L140: i__1 = *liw - 1; ma27q_(n, &a[1], la, &iw[2], &i__1, &w[1], maxfrt, &rhs[1], &iw1[1], & nblk, &latop, &icntl[1]); /* BACK SUBSTITUTION. */ i__1 = *liw - 1; ma27r_(n, &a[1], la, &iw[2], &i__1, &w[1], maxfrt, &rhs[1], &iw1[1], & nblk, &latop, &icntl[1]); L150: if (icntl[3] <= 0 || icntl[2] <= 0) { goto L170; } /* PRINT OUTPUT PARAMETERS. */ io___81.ciunit = icntl[2]; s_wsfe(&io___81); e_wsfe(); if (*n > 0) { io___82.ciunit = icntl[2]; s_wsfe(&io___82); i__1 = k; for (i__ = 1; i__ <= i__1; ++i__) { do_fio(&c__1, (char *)&rhs[i__], (ftnlen)sizeof(real)); } e_wsfe(); } L170: return 0; } /* ma27c_ */ /* Subroutine */ int ma27g_(integer *n, integer *nz, integer *irn, integer * icn, integer *iw, integer *lw, integer *ipe, integer *iq, integer * flag__, integer *iwfr, integer *icntl, integer *info) { /* Format strings */ static char fmt_60[] = "(\002 *** WARNING MESSAGE FROM SUBROUTINE MA27" "A\002,\002 *** INFO(1) =\002,i2)"; static char fmt_70[] = "(i6,\002TH NON-ZERO (IN ROW\002,i6,\002 AND COLU" "MN\002,i6,\002) IGNORED\002)"; /* System generated locals */ integer i__1, i__2; /* Builtin functions */ integer s_wsfe(cilist *), do_fio(integer *, char *, ftnlen), e_wsfe(void); /* Local variables */ static integer i__, j, k, l, k1, k2, n1, id, jn, lr, last, ndup; /* Fortran I/O blocks */ static cilist io___87 = { 0, 0, 0, fmt_60, 0 }; static cilist io___88 = { 0, 0, 0, fmt_70, 0 }; /* SORT PRIOR TO CALLING ANALYSIS ROUTINE MA27H/HD. */ /* GIVEN THE POSITIONS OF THE OFF-DIAGONAL NON-ZEROS OF A SYMMETRIC */ /* MATRIX, CONSTRUCT THE SPARSITY PATTERN OF THE OFF-DIAGONAL */ /* PART OF THE WHOLE MATRIX (UPPER AND LOWER TRIANGULAR PARTS). */ /* EITHER ONE OF A PAIR (I,J),(J,I) MAY BE USED TO REPRESENT */ /* THE PAIR. DIAGONAL ELEMENTS AND DUPLICATES ARE IGNORED. */ /* N MUST BE SET TO THE MATRIX ORDER. IT IS NOT ALTERED. */ /* NZ MUST BE SET TO THE NUMBER OF NON-ZEROS INPUT. IT IS NOT */ /* ALTERED. */ /* IRN(I),I=1,2,...,NZ MUST BE SET TO THE ROW NUMBERS OF THE */ /* NON-ZEROS ON INPUT. IT IS NOT ALTERED UNLESS IT IS EQUIVALENCED */ /* TO IW (SEE DESCRIPTION OF IW). */ /* ICN(I),I=1,2,...,NZ MUST BE SET TO THE COLUMN NUMBERS OF THE */ /* NON-ZEROS ON INPUT. IT IS NOT ALTERED UNLESS IT IS EQUIVALENCED */ /* TO IW (SEE DESCRIPTION OF IW). */ /* IW NEED NOT BE SET ON INPUT. ON OUTPUT IT CONTAINS LISTS OF */ /* COLUMN INDICES, EACH LIST BEING HEADED BY ITS LENGTH. */ /* IRN(1) MAY BE EQUIVALENCED TO IW(1) AND ICN(1) MAY BE */ /* EQUIVALENCED TO IW(K), WHERE K.GT.NZ. */ /* LW MUST BE SET TO THE LENGTH OF IW. IT MUST BE AT LEAST 2*NZ+N. */ /* IT IS NOT ALTERED. */ /* IPE NEED NOT BE SET ON INPUT. ON OUTPUT IPE(I) POINTS TO THE START OF */ /* THE ENTRY IN IW FOR ROW I, OR IS ZERO IF THERE IS NO ENTRY. */ /* IQ NEED NOT BE SET. ON OUTPUT IQ(I),I=1,N CONTAINS THE NUMBER OF */ /* OFF-DIAGONAL N0N-ZEROS IN ROW I INCLUDING DUPLICATES. */ /* FLAG IS USED FOR WORKSPACE TO HOLD FLAGS TO PERMIT DUPLICATE ENTRIES */ /* TO BE IDENTIFIED QUICKLY. */ /* IWFR NEED NOT BE SET ON INPUT. ON OUTPUT IT POINTS TO THE FIRST */ /* UNUSED LOCATION IN IW. */ /* ICNTL is an INTEGER array of length 30, see MA27A/AD. */ /* INFO is an INTEGER array of length 20, see MA27A/AD. */ /* .. Scalar Arguments .. */ /* .. */ /* .. Array Arguments .. */ /* .. */ /* .. Local Scalars .. */ /* .. */ /* .. Executable Statements .. */ /* INITIALIZE INFO(2) AND COUNT IN IPE THE */ /* NUMBERS OF NON-ZEROS IN THE ROWS AND MOVE ROW AND COLUMN */ /* NUMBERS INTO IW. */ /* Parameter adjustments */ --flag__; --iq; --ipe; --irn; --icn; --iw; --icntl; --info; /* Function Body */ info[2] = 0; i__1 = *n; for (i__ = 1; i__ <= i__1; ++i__) { ipe[i__] = 0; /* L10: */ } lr = *nz; if (*nz == 0) { goto L120; } i__1 = *nz; for (k = 1; k <= i__1; ++k) { i__ = irn[k]; j = icn[k]; if (i__ < j) { if (i__ >= 1 && j <= *n) { goto L90; } } else if (i__ > j) { if (j >= 1 && i__ <= *n) { goto L90; } } else { if (i__ >= 1 && i__ <= *n) { goto L80; } } ++info[2]; info[1] = 1; if (info[2] <= 1 && icntl[2] > 0) { io___87.ciunit = icntl[2]; s_wsfe(&io___87); do_fio(&c__1, (char *)&info[1], (ftnlen)sizeof(integer)); e_wsfe(); } if (info[2] <= 10 && icntl[2] > 0) { io___88.ciunit = icntl[2]; s_wsfe(&io___88); do_fio(&c__1, (char *)&k, (ftnlen)sizeof(integer)); do_fio(&c__1, (char *)&i__, (ftnlen)sizeof(integer)); do_fio(&c__1, (char *)&j, (ftnlen)sizeof(integer)); e_wsfe(); } L80: i__ = 0; j = 0; goto L100; L90: ++ipe[i__]; ++ipe[j]; L100: iw[k] = j; ++lr; iw[lr] = i__; /* L110: */ } /* ACCUMULATE ROW COUNTS TO GET POINTERS TO ROW STARTS IN BOTH IPE AND IQ */ /* AND INITIALIZE FLAG */ L120: iq[1] = 1; n1 = *n - 1; if (n1 <= 0) { goto L140; } i__1 = n1; for (i__ = 1; i__ <= i__1; ++i__) { flag__[i__] = 0; if (ipe[i__] == 0) { ipe[i__] = -1; } iq[i__ + 1] = ipe[i__] + iq[i__] + 1; ipe[i__] = iq[i__]; /* L130: */ } L140: last = ipe[*n] + iq[*n]; flag__[*n] = 0; if (lr >= last) { goto L160; } k1 = lr + 1; i__1 = last; for (k = k1; k <= i__1; ++k) { iw[k] = 0; /* L150: */ } L160: ipe[*n] = iq[*n]; *iwfr = last + 1; if (*nz == 0) { goto L230; } /* RUN THROUGH PUTTING THE MATRIX ELEMENTS IN THE RIGHT PLACE */ /* BUT WITH SIGNS INVERTED. IQ IS USED FOR HOLDING RUNNING POINTERS */ /* AND IS LEFT HOLDING POINTERS TO ROW ENDS. */ i__1 = *nz; for (k = 1; k <= i__1; ++k) { j = iw[k]; if (j <= 0) { goto L220; } l = k; iw[k] = 0; i__2 = *nz; for (id = 1; id <= i__2; ++id) { if (l > *nz) { goto L170; } l += *nz; goto L180; L170: l -= *nz; L180: i__ = iw[l]; iw[l] = 0; if (i__ < j) { goto L190; } l = iq[j] + 1; iq[j] = l; jn = iw[l]; iw[l] = -i__; goto L200; L190: l = iq[i__] + 1; iq[i__] = l; jn = iw[l]; iw[l] = -j; L200: j = jn; if (j <= 0) { goto L220; } /* L210: */ } L220: ; } /* RUN THROUGH RESTORING SIGNS, REMOVING DUPLICATES AND SETTING THE */ /* MATE OF EACH NON-ZERO. */ /* NDUP COUNTS THE NUMBER OF DUPLICATE ELEMENTS. */ L230: ndup = 0; i__1 = *n; for (i__ = 1; i__ <= i__1; ++i__) { k1 = ipe[i__] + 1; k2 = iq[i__]; if (k1 <= k2) { goto L240; } /* ROW IS EMPTY. SET POINTER TO ZERO. */ ipe[i__] = 0; iq[i__] = 0; goto L280; /* ON ENTRY TO THIS LOOP FLAG(J).LT.I FOR J=1,2,...,N. DURING THE LOOP */ /* FLAG(J) IS SET TO I IF A NON-ZERO IN COLUMN J IS FOUND. THIS */ /* PERMITS DUPLICATES TO BE RECOGNIZED QUICKLY. */ L240: i__2 = k2; for (k = k1; k <= i__2; ++k) { j = -iw[k]; if (j <= 0) { goto L270; } l = iq[j] + 1; iq[j] = l; iw[l] = i__; iw[k] = j; if (flag__[j] != i__) { goto L250; } ++ndup; iw[l] = 0; iw[k] = 0; L250: flag__[j] = i__; /* L260: */ } L270: iq[i__] -= ipe[i__]; if (ndup == 0) { iw[k1 - 1] = iq[i__]; } L280: ; } if (ndup == 0) { goto L310; } /* COMPRESS IW TO REMOVE DUMMY ENTRIES CAUSED BY DUPLICATES. */ *iwfr = 1; i__1 = *n; for (i__ = 1; i__ <= i__1; ++i__) { k1 = ipe[i__] + 1; if (k1 == 1) { goto L300; } k2 = iq[i__] + ipe[i__]; l = *iwfr; ipe[i__] = *iwfr; ++(*iwfr); i__2 = k2; for (k = k1; k <= i__2; ++k) { if (iw[k] == 0) { goto L290; } iw[*iwfr] = iw[k]; ++(*iwfr); L290: ; } iw[l] = *iwfr - l - 1; L300: ; } L310: return 0; } /* ma27g_ */ /* Subroutine */ int ma27h_(integer *n, integer *ipe, integer *iw, integer * lw, integer *iwfr, integer *nv, integer *nxt, integer *lst, integer * ipd, integer *flag__, integer *iovflo, integer *ncmpa, real *fratio) { /* System generated locals */ integer i__1, i__2, i__3, i__4, i__5; real r__1; /* Builtin functions */ integer i_nint(real *); /* Local variables */ static integer i__, k, l, k1, k2, id, ie, ke, md, me, ip, jp, kp, is, js, ks, ln, ls, ml, ms, np, ns, jp1, jp2, kp0, kp1, kp2, np0, idl, idn, len, nel, nflg; extern /* Subroutine */ int ma27u_(integer *, integer *, integer *, integer *, integer *, integer *); static integer lwfr, root, limit, nvpiv, nvroot; /* ANALYSIS SUBROUTINE */ /* GIVEN REPRESENTATION OF THE WHOLE MATRIX (EXCLUDING DIAGONAL) */ /* PERFORM MINIMUM DEGREE ORDERING, CONSTRUCTING TREE POINTERS. */ /* IT WORKS WITH SUPERVARIABLES WHICH ARE COLLECTIONS OF ONE OR MORE */ /* VARIABLES, STARTING WITH SUPERVARIABLE I CONTAINING VARIABLE I FOR */ /* I = 1,2,...,N. ALL VARIABLES IN A SUPERVARIABLE ARE ELIMINATED */ /* TOGETHER. EACH SUPERVARIABLE HAS AS NUMERICAL NAME THAT OF ONE */ /* OF ITS VARIABLES (ITS PRINCIPAL VARIABLE). */ /* N MUST BE SET TO THE MATRIX ORDER. IT IS NOT ALTERED. */ /* IPE(I) MUST BE SET TO POINT TO THE POSITION IN IW OF THE */ /* START OF ROW I OR HAVE THE VALUE ZERO IF ROW I HAS NO OFF- */ /* DIAGONAL NON-ZEROS. DURING EXECUTION IT IS USED AS FOLLOWS. IF */ /* SUPERVARIABLE I IS ABSORBED INTO SUPERVARIABLE J THEN IPE(I)=-J. */ /* IF SUPERVARIABLE I IS ELIMINATED THEN IPE(I) EITHER POINTS TO THE */ /* LIST OF SUPERVARIABLES FOR CREATED ELEMENT I OR IS ZERO IF */ /* THE CREATED ELEMENT IS NULL. IF ELEMENT I */ /* IS ABSORBED INTO ELEMENT J THEN IPE(I)=-J. */ /* IW MUST BE SET ON ENTRY TO HOLD LISTS OF VARIABLES BY */ /* ROWS, EACH LIST BEING HEADED BY ITS LENGTH. */ /* DURING EXECUTION THESE LISTS ARE REVISED AND HOLD */ /* LISTS OF ELEMENTS AND SUPERVARIABLES. THE ELEMENTS */ /* ALWAYS HEAD THE LISTS. WHEN A SUPERVARIABLE */ /* IS ELIMINATED ITS LIST IS REPLACED BY A LIST OF SUPERVARIABLES */ /* IN THE NEW ELEMENT. */ /* LW MUST BE SET TO THE LENGTH OF IW. IT IS NOT ALTERED. */ /* IWFR MUST BE SET TO THE POSITION IN IW OF THE FIRST FREE VARIABLE. */ /* IT IS REVISED DURING EXECUTION AND CONTINUES TO HAVE THIS MEANING. */ /* NV(JS) NEED NOT BE SET. DURING EXECUTION IT IS ZERO IF */ /* JS IS NOT A PRINCIPAL VARIABLE AND IF IT IS IT HOLDS */ /* THE NUMBER OF VARIABLES IN SUPERVARIABLE JS. FOR ELIMINATED */ /* VARIABLES IT IS SET TO THE DEGREE AT THE TIME OF ELIMINATION. */ /* NXT(JS) NEED NOT BE SET. DURING EXECUTION IT IS THE NEXT */ /* SUPERVARIABLE HAVING THE SAME DEGREE AS JS, OR ZERO */ /* IF IT IS LAST IN ITS LIST. */ /* LST(JS) NEED NOT BE SET. DURING EXECUTION IT IS THE */ /* LAST SUPERVARIABLE HAVING THE SAME DEGREE AS JS OR */ /* -(ITS DEGREE) IF IT IS FIRST IN ITS LIST. */ /* IPD(ID) NEED NOT BE SET. DURING EXECUTION IT */ /* IS THE FIRST SUPERVARIABLE WITH DEGREE ID OR ZERO */ /* IF THERE ARE NONE. */ /* FLAG IS USED AS WORKSPACE FOR ELEMENT AND SUPERVARIABLE FLAGS. */ /* WHILE THE CODE IS FINDING THE DEGREE OF SUPERVARIABLE IS */ /* FLAG HAS THE FOLLOWING VALUES. */ /* A) FOR THE CURRENT PIVOT/NEW ELEMENT ME */ /* FLAG(ME)=-1 */ /* B) FOR VARIABLES JS */ /* FLAG(JS)=-1 IF JS IS NOT A PRINCIPAL VARIABLE */ /* FLAG(JS)=0 IF JS IS A SUPERVARIABLE IN THE NEW ELEMENT */ /* FLAG(JS)=NFLG IF JS IS A SUPERVARIABLE NOT IN THE NEW */ /* ELEMENT THAT HAS BEEN COUNTED IN THE DEGREE */ /* CALCULATION */ /* FLAG(JS).GT.NFLG IF JS IS A SUPERVARIABLE NOT IN THE NEW */ /* ELEMENT THAT HAS NOT BEEN COUNTED IN THE DEGREE */ /* CALCULATION */ /* C) FOR ELEMENTS IE */ /* FLAG(IE)=-1 IF ELEMENT IE HAS BEEN MERGED INTO ANOTHER */ /* FLAG(IE)=-NFLG IF ELEMENT IE HAS BEEN USED IN THE DEGREE */ /* CALCULATION FOR IS. */ /* FLAG(IE).LT.-NFLG IF ELEMENT IE HAS NOT BEEN USED IN THE */ /* DEGREE CALCULATION FOR IS */ /* IOVFLO see ICNTL(4) in MA27A/AD. */ /* NCMPA see INFO(11) in MA27A/AD. */ /* FRATIO see CNTL(2) in MA27A/AD. */ /* .. Scalar Arguments .. */ /* .. */ /* .. Array Arguments .. */ /* .. */ /* .. Local Scalars .. */ /* LIMIT Limit on number of variables for putting node in root. */ /* NVROOT Number of variables in the root node */ /* ROOT Index of the root node (N+1 if none chosen yet). */ /* .. */ /* .. External Subroutines .. */ /* .. */ /* .. Intrinsic Functions .. */ /* .. */ /* If a column of the reduced matrix has relative density greater than */ /* CNTL(2), it is forced into the root. All such columns are taken to */ /* have sparsity pattern equal to their merged patterns, so the fill */ /* and operation counts may be overestimated. */ /* IS,JS,KS,LS,MS,NS ARE USED TO REFER TO SUPERVARIABLES. */ /* IE,JE,KE ARE USED TO REFER TO ELEMENTS. */ /* IP,JP,KP,K,NP ARE USED TO POINT TO LISTS OF ELEMENTS. */ /* OR SUPERVARIABLES. */ /* ID IS USED FOR THE DEGREE OF A SUPERVARIABLE. */ /* MD IS USED FOR THE CURRENT MINIMUM DEGREE. */ /* IDN IS USED FOR THE NO. OF VARIABLES IN A NEWLY CREATED ELEMENT */ /* NEL IS USED TO HOLD THE NO. OF VARIABLES THAT HAVE BEEN */ /* ELIMINATED. */ /* ME=MS IS THE NAME OF THE SUPERVARIABLE ELIMINATED AND */ /* OF THE ELEMENT CREATED IN THE MAIN LOOP. */ /* NFLG IS USED FOR THE CURRENT FLAG VALUE IN ARRAY FLAG. IT STARTS */ /* WITH THE VALUE IOVFLO AND IS REDUCED BY 1 EACH TIME IT IS USED */ /* UNTIL IT HAS THE VALUE 2 WHEN IT IS RESET TO THE VALUE IOVFLO. */ /* .. Executable Statements .. */ /* INITIALIZATIONS */ /* Parameter adjustments */ --flag__; --ipd; --lst; --nxt; --nv; --ipe; --iw; /* Function Body */ i__1 = *n; for (i__ = 1; i__ <= i__1; ++i__) { ipd[i__] = 0; nv[i__] = 1; flag__[i__] = *iovflo; /* L10: */ } md = 1; *ncmpa = 0; nflg = *iovflo; nel = 0; root = *n + 1; nvroot = 0; /* LINK TOGETHER VARIABLES HAVING SAME DEGREE */ i__1 = *n; for (is = 1; is <= i__1; ++is) { k = ipe[is]; if (k <= 0) { goto L20; } id = iw[k] + 1; ns = ipd[id]; if (ns > 0) { lst[ns] = is; } nxt[is] = ns; ipd[id] = is; lst[is] = -id; goto L30; /* WE HAVE A VARIABLE THAT CAN BE ELIMINATED AT ONCE BECAUSE THERE IS */ /* NO OFF-DIAGONAL NON-ZERO IN ITS ROW. */ L20: ++nel; flag__[is] = -1; nxt[is] = 0; lst[is] = 0; L30: ; } /* START OF MAIN LOOP */ i__1 = *n; for (ml = 1; ml <= i__1; ++ml) { /* LEAVE LOOP IF ALL VARIABLES HAVE BEEN ELIMINATED. */ if (nel + nvroot + 1 >= *n) { goto L350; } /* FIND NEXT SUPERVARIABLE FOR ELIMINATION. */ i__2 = *n; for (id = md; id <= i__2; ++id) { ms = ipd[id]; if (ms > 0) { goto L50; } /* L40: */ } L50: md = id; /* NVPIV HOLDS THE NUMBER OF VARIABLES IN THE PIVOT. */ nvpiv = nv[ms]; /* REMOVE CHOSEN VARIABLE FROM LINKED LIST */ ns = nxt[ms]; nxt[ms] = 0; lst[ms] = 0; if (ns > 0) { lst[ns] = -id; } ipd[id] = ns; me = ms; nel += nvpiv; /* IDN HOLDS THE DEGREE OF THE NEW ELEMENT. */ idn = 0; /* RUN THROUGH THE LIST OF THE PIVOTAL SUPERVARIABLE, SETTING TREE */ /* POINTERS AND CONSTRUCTING NEW LIST OF SUPERVARIABLES. */ /* KP IS A POINTER TO THE CURRENT POSITION IN THE OLD LIST. */ kp = ipe[me]; flag__[ms] = -1; /* IP POINTS TO THE START OF THE NEW LIST. */ ip = *iwfr; /* LEN HOLDS THE LENGTH OF THE LIST ASSOCIATED WITH THE PIVOT. */ len = iw[kp]; i__2 = len; for (kp1 = 1; kp1 <= i__2; ++kp1) { ++kp; ke = iw[kp]; /* JUMP IF KE IS AN ELEMENT THAT HAS NOT BEEN MERGED INTO ANOTHER. */ if (flag__[ke] <= -2) { goto L60; } /* JUMP IF KE IS AN ELEMENT THAT HAS BEEN MERGED INTO ANOTHER OR IS */ /* A SUPERVARIABLE THAT HAS BEEN ELIMINATED. */ if (flag__[ke] <= 0) { if (ipe[ke] != -root) { goto L140; } /* KE has been merged into the root */ ke = root; if (flag__[ke] <= 0) { goto L140; } } /* WE HAVE A SUPERVARIABLE. PREPARE TO SEARCH REST OF LIST. */ jp = kp - 1; ln = len - kp1 + 1; ie = ms; goto L70; /* SEARCH VARIABLE LIST OF ELEMENT KE, USING JP AS A POINTER TO IT. */ L60: ie = ke; jp = ipe[ie]; ln = iw[jp]; /* SEARCH FOR DIFFERENT SUPERVARIABLES AND ADD THEM TO THE NEW LIST, */ /* COMPRESSING WHEN NECESSARY. THIS LOOP IS EXECUTED ONCE FOR */ /* EACH ELEMENT IN THE LIST AND ONCE FOR ALL THE SUPERVARIABLES */ /* IN THE LIST. */ L70: i__3 = ln; for (jp1 = 1; jp1 <= i__3; ++jp1) { ++jp; is = iw[jp]; /* JUMP IF IS IS NOT A PRINCIPAL VARIABLE OR HAS ALREADY BEEN COUNTED. */ if (flag__[is] <= 0) { if (ipe[is] == -root) { /* IS has been merged into the root */ is = root; iw[jp] = root; if (flag__[is] <= 0) { goto L130; } } else { goto L130; } } flag__[is] = 0; if (*iwfr < *lw) { goto L100; } /* PREPARE FOR COMPRESSING IW BY ADJUSTING POINTERS AND */ /* LENGTHS SO THAT THE LISTS BEING SEARCHED IN THE INNER AND OUTER */ /* LOOPS CONTAIN ONLY THE REMAINING ENTRIES. */ ipe[ms] = kp; iw[kp] = len - kp1; ipe[ie] = jp; iw[jp] = ln - jp1; /* COMPRESS IW */ i__4 = ip - 1; ma27u_(n, &ipe[1], &iw[1], &i__4, &lwfr, ncmpa); /* COPY NEW LIST FORWARD */ jp2 = *iwfr - 1; *iwfr = lwfr; if (ip > jp2) { goto L90; } i__4 = jp2; for (jp = ip; jp <= i__4; ++jp) { iw[*iwfr] = iw[jp]; ++(*iwfr); /* L80: */ } /* ADJUST POINTERS FOR THE NEW LIST AND THE LISTS BEING SEARCHED. */ L90: ip = lwfr; jp = ipe[ie]; kp = ipe[me]; /* STORE IS IN NEW LIST. */ L100: iw[*iwfr] = is; idn += nv[is]; ++(*iwfr); /* REMOVE IS FROM DEGREE LINKED LIST */ ls = lst[is]; lst[is] = 0; ns = nxt[is]; nxt[is] = 0; if (ns > 0) { lst[ns] = ls; } if (ls < 0) { ls = -ls; ipd[ls] = ns; } else if (ls > 0) { nxt[ls] = ns; } L130: ; } /* JUMP IF WE HAVE JUST BEEN SEARCHING THE VARIABLES AT THE END OF */ /* THE LIST OF THE PIVOT. */ if (ie == ms) { goto L150; } /* SET TREE POINTER AND FLAG TO INDICATE ELEMENT IE IS ABSORBED INTO */ /* NEW ELEMENT ME. */ ipe[ie] = -me; flag__[ie] = -1; L140: ; } /* STORE THE DEGREE OF THE PIVOT. */ L150: nv[ms] = idn + nvpiv; /* JUMP IF NEW ELEMENT IS NULL. */ if (*iwfr == ip) { goto L330; } k1 = ip; k2 = *iwfr - 1; /* RUN THROUGH NEW LIST OF SUPERVARIABLES REVISING EACH ASSOCIATED LIST, */ /* RECALCULATING DEGREES AND REMOVING DUPLICATES. */ r__1 = *fratio * (*n - nel); limit = i_nint(&r__1); i__2 = k2; for (k = k1; k <= i__2; ++k) { is = iw[k]; if (is == root) { goto L310; } if (nflg > 2) { goto L170; } /* RESET FLAG VALUES TO +/-IOVFLO. */ i__3 = *n; for (i__ = 1; i__ <= i__3; ++i__) { if (flag__[i__] > 0) { flag__[i__] = *iovflo; } if (flag__[i__] <= -2) { flag__[i__] = -(*iovflo); } /* L160: */ } nflg = *iovflo; /* REDUCE NFLG BY ONE TO CATER FOR THIS SUPERVARIABLE. */ L170: --nflg; /* BEGIN WITH THE DEGREE OF THE NEW ELEMENT. ITS VARIABLES MUST ALWAYS */ /* BE COUNTED DURING THE DEGREE CALCULATION AND THEY ARE ALREADY */ /* FLAGGED WITH THE VALUE 0. */ id = idn; /* RUN THROUGH THE LIST ASSOCIATED WITH SUPERVARIABLE IS */ kp1 = ipe[is] + 1; /* NP POINTS TO THE NEXT ENTRY IN THE REVISED LIST. */ np = kp1; kp2 = iw[kp1 - 1] + kp1 - 1; i__3 = kp2; for (kp = kp1; kp <= i__3; ++kp) { ke = iw[kp]; /* TEST WHETHER KE IS AN ELEMENT, A REDUNDANT ENTRY OR A SUPERVARIABLE. */ if (flag__[ke] == -1) { if (ipe[ke] != -root) { goto L220; } /* KE has been merged into the root */ ke = root; iw[kp] = root; if (flag__[ke] == -1) { goto L220; } } if (flag__[ke] >= 0) { goto L230; } /* SEARCH LIST OF ELEMENT KE, REVISING THE DEGREE WHEN NEW VARIABLES */ /* FOUND. */ jp1 = ipe[ke] + 1; jp2 = iw[jp1 - 1] + jp1 - 1; idl = id; i__4 = jp2; for (jp = jp1; jp <= i__4; ++jp) { js = iw[jp]; /* JUMP IF JS HAS ALREADY BEEN COUNTED. */ if (flag__[js] <= nflg) { goto L190; } id += nv[js]; flag__[js] = nflg; L190: ; } /* JUMP IF ONE OR MORE NEW SUPERVARIABLES WERE FOUND. */ if (id > idl) { goto L210; } /* CHECK WHETHER EVERY VARIABLE OF ELEMENT KE IS IN NEW ELEMENT ME. */ i__4 = jp2; for (jp = jp1; jp <= i__4; ++jp) { js = iw[jp]; if (flag__[js] != 0) { goto L210; } /* L200: */ } /* SET TREE POINTER AND FLAG TO INDICATE THAT ELEMENT KE IS ABSORBED */ /* INTO NEW ELEMENT ME. */ ipe[ke] = -me; flag__[ke] = -1; goto L220; /* STORE ELEMENT KE IN THE REVISED LIST FOR SUPERVARIABLE IS AND FLAG IT. */ L210: iw[np] = ke; flag__[ke] = -nflg; ++np; L220: ; } np0 = np; goto L250; /* TREAT THE REST OF THE LIST ASSOCIATED WITH SUPERVARIABLE IS. IT */ /* CONSISTS ENTIRELY OF SUPERVARIABLES. */ L230: kp0 = kp; np0 = np; i__3 = kp2; for (kp = kp0; kp <= i__3; ++kp) { ks = iw[kp]; if (flag__[ks] <= nflg) { if (ipe[ks] == -root) { ks = root; iw[kp] = root; if (flag__[ks] <= nflg) { goto L240; } } else { goto L240; } } /* ADD TO DEGREE, FLAG SUPERVARIABLE KS AND ADD IT TO NEW LIST. */ id += nv[ks]; flag__[ks] = nflg; iw[np] = ks; ++np; L240: ; } /* MOVE FIRST SUPERVARIABLE TO END OF LIST, MOVE FIRST ELEMENT TO END */ /* OF ELEMENT PART OF LIST AND ADD NEW ELEMENT TO FRONT OF LIST. */ L250: if (id >= limit) { goto L295; } iw[np] = iw[np0]; iw[np0] = iw[kp1]; iw[kp1] = me; /* STORE THE NEW LENGTH OF THE LIST. */ iw[kp1 - 1] = np - kp1 + 1; /* CHECK WHETHER ROW IS IS IDENTICAL TO ANOTHER BY LOOKING IN LINKED */ /* LIST OF SUPERVARIABLES WITH DEGREE ID AT THOSE WHOSE LISTS HAVE */ /* FIRST ENTRY ME. NOTE THAT THOSE CONTAINING ME COME FIRST SO THE */ /* SEARCH CAN BE TERMINATED WHEN A LIST NOT STARTING WITH ME IS */ /* FOUND. */ js = ipd[id]; i__3 = *n; for (l = 1; l <= i__3; ++l) { if (js <= 0) { goto L300; } kp1 = ipe[js] + 1; if (iw[kp1] != me) { goto L300; } /* JS HAS SAME DEGREE AND IS ACTIVE. CHECK IF IDENTICAL TO IS. */ kp2 = kp1 - 1 + iw[kp1 - 1]; i__4 = kp2; for (kp = kp1; kp <= i__4; ++kp) { ie = iw[kp]; /* JUMP IF IE IS A SUPERVARIABLE OR AN ELEMENT NOT IN THE LIST OF IS. */ if ((i__5 = flag__[ie], abs(i__5)) > nflg) { goto L270; } /* L260: */ } goto L290; L270: js = nxt[js]; /* L280: */ } /* SUPERVARIABLE AMALGAMATION. ROW IS IS IDENTICAL TO ROW JS. */ /* REGARD ALL VARIABLES IN THE TWO SUPERVARIABLES AS BEING IN IS. SET */ /* TREE POINTER, FLAG AND NV ENTRIES. */ L290: ipe[js] = -is; nv[is] += nv[js]; nv[js] = 0; flag__[js] = -1; /* REPLACE JS BY IS IN LINKED LIST. */ ns = nxt[js]; ls = lst[js]; if (ns > 0) { lst[ns] = is; } if (ls > 0) { nxt[ls] = is; } lst[is] = ls; nxt[is] = ns; lst[js] = 0; nxt[js] = 0; if (ipd[id] == js) { ipd[id] = is; } goto L310; /* Treat IS as full. Merge it into the root node. */ L295: if (nvroot == 0) { root = is; ipe[is] = 0; } else { iw[k] = root; ipe[is] = -root; nv[root] += nv[is]; nv[is] = 0; flag__[is] = -1; } nvroot = nv[root]; goto L310; /* INSERT IS INTO LINKED LIST OF SUPERVARIABLES OF SAME DEGREE. */ L300: ns = ipd[id]; if (ns > 0) { lst[ns] = is; } nxt[is] = ns; ipd[id] = is; lst[is] = -id; md = min(md,id); L310: ; } /* RESET FLAGS FOR SUPERVARIABLES IN NEWLY CREATED ELEMENT AND */ /* REMOVE THOSE ABSORBED INTO OTHERS. */ i__2 = k2; for (k = k1; k <= i__2; ++k) { is = iw[k]; if (nv[is] == 0) { goto L320; } flag__[is] = nflg; iw[ip] = is; ++ip; L320: ; } *iwfr = k1; flag__[me] = -nflg; /* MOVE FIRST ENTRY TO END TO MAKE ROOM FOR LENGTH. */ iw[ip] = iw[k1]; iw[k1] = ip - k1; /* SET POINTER FOR NEW ELEMENT AND RESET IWFR. */ ipe[me] = k1; *iwfr = ip + 1; goto L335; L330: ipe[me] = 0; L335: /* L340: */ ; } /* Absorb any remaining variables into the root */ L350: i__1 = *n; for (is = 1; is <= i__1; ++is) { if (nxt[is] != 0 || lst[is] != 0) { if (nvroot == 0) { root = is; ipe[is] = 0; } else { ipe[is] = -root; } nvroot += nv[is]; nv[is] = 0; } /* L360: */ } /* Link any remaining elements to the root */ i__1 = *n; for (ie = 1; ie <= i__1; ++ie) { if (ipe[ie] > 0) { ipe[ie] = -root; } /* L370: */ } if (nvroot > 0) { nv[root] = nvroot; } return 0; } /* ma27h_ */ /* Subroutine */ int ma27u_(integer *n, integer *ipe, integer *iw, integer * lw, integer *iwfr, integer *ncmpa) { /* System generated locals */ integer i__1, i__2; /* Local variables */ static integer i__, k, k1, k2, ir, lwfr; /* COMPRESS LISTS HELD BY MA27H/HD AND MA27K/KD IN IW AND ADJUST POINTERS */ /* IN IPE TO CORRESPOND. */ /* N IS THE MATRIX ORDER. IT IS NOT ALTERED. */ /* IPE(I) POINTS TO THE POSITION IN IW OF THE START OF LIST I OR IS */ /* ZERO IF THERE IS NO LIST I. ON EXIT IT POINTS TO THE NEW POSITION. */ /* IW HOLDS THE LISTS, EACH HEADED BY ITS LENGTH. ON OUTPUT THE SAME */ /* LISTS ARE HELD, BUT THEY ARE NOW COMPRESSED TOGETHER. */ /* LW HOLDS THE LENGTH OF IW. IT IS NOT ALTERED. */ /* IWFR NEED NOT BE SET ON ENTRY. ON EXIT IT POINTS TO THE FIRST FREE */ /* LOCATION IN IW. */ /* ON RETURN IT IS SET TO THE FIRST FREE LOCATION IN IW. */ /* NCMPA see INFO(11) in MA27A/AD. */ /* .. Scalar Arguments .. */ /* .. */ /* .. Array Arguments .. */ /* .. */ /* .. Local Scalars .. */ /* .. */ /* .. Executable Statements .. */ /* Parameter adjustments */ --ipe; --iw; /* Function Body */ ++(*ncmpa); /* PREPARE FOR COMPRESSING BY STORING THE LENGTHS OF THE */ /* LISTS IN IPE AND SETTING THE FIRST ENTRY OF EACH LIST TO */ /* -(LIST NUMBER). */ i__1 = *n; for (i__ = 1; i__ <= i__1; ++i__) { k1 = ipe[i__]; if (k1 <= 0) { goto L10; } ipe[i__] = iw[k1]; iw[k1] = -i__; L10: ; } /* COMPRESS */ /* IWFR POINTS JUST BEYOND THE END OF THE COMPRESSED FILE. */ /* LWFR POINTS JUST BEYOND THE END OF THE UNCOMPRESSED FILE. */ *iwfr = 1; lwfr = *iwfr; i__1 = *n; for (ir = 1; ir <= i__1; ++ir) { if (lwfr > *lw) { goto L70; } /* SEARCH FOR THE NEXT NEGATIVE ENTRY. */ i__2 = *lw; for (k = lwfr; k <= i__2; ++k) { if (iw[k] < 0) { goto L30; } /* L20: */ } goto L70; /* PICK UP ENTRY NUMBER, STORE LENGTH IN NEW POSITION, SET NEW POINTER */ /* AND PREPARE TO COPY LIST. */ L30: i__ = -iw[k]; iw[*iwfr] = ipe[i__]; ipe[i__] = *iwfr; k1 = k + 1; k2 = k + iw[*iwfr]; ++(*iwfr); if (k1 > k2) { goto L50; } /* COPY LIST TO NEW POSITION. */ i__2 = k2; for (k = k1; k <= i__2; ++k) { iw[*iwfr] = iw[k]; ++(*iwfr); /* L40: */ } L50: lwfr = k2 + 1; /* L60: */ } L70: return 0; } /* ma27u_ */ /* Subroutine */ int ma27j_(integer *n, integer *nz, integer *irn, integer * icn, integer *perm, integer *iw, integer *lw, integer *ipe, integer * iq, integer *flag__, integer *iwfr, integer *icntl, integer *info) { /* Format strings */ static char fmt_60[] = "(\002 *** WARNING MESSAGE FROM SUBROUTINE MA27" "A\002,\002 *** INFO(1) =\002,i2)"; static char fmt_70[] = "(i6,\002TH NON-ZERO (IN ROW\002,i6,\002 AND COLU" "MN\002,i6,\002) IGNORED\002)"; /* System generated locals */ integer i__1, i__2; /* Builtin functions */ integer s_wsfe(cilist *), do_fio(integer *, char *, ftnlen), e_wsfe(void); /* Local variables */ static integer i__, j, k, l, k1, k2, id, in, len, lbig, jdummy; /* Fortran I/O blocks */ static cilist io___144 = { 0, 0, 0, fmt_60, 0 }; static cilist io___145 = { 0, 0, 0, fmt_70, 0 }; /* SORT PRIOR TO CALLING ANALYSIS ROUTINE MA27K/KD. */ /* GIVEN THE POSITIONS OF THE OFF-DIAGONAL NON-ZEROS OF A SYMMETRIC */ /* MATRIX AND A PERMUTATION, CONSTRUCT THE SPARSITY PATTERN */ /* OF THE STRICTLY UPPER TRIANGULAR PART OF THE PERMUTED MATRIX. */ /* EITHER ONE OF A PAIR (I,J),(J,I) MAY BE USED TO REPRESENT */ /* THE PAIR. DIAGONAL ELEMENTS ARE IGNORED. NO CHECK IS MADE */ /* FOR DUPLICATE ELEMENTS UNLESS ANY ROW HAS MORE THAN ICNTL(4) */ /* NON-ZEROS, IN WHICH CASE DUPLICATES ARE REMOVED. */ /* N MUST BE SET TO THE MATRIX ORDER. IT IS NOT ALTERED. */ /* NZ MUST BE SET TO THE NUMBER OF NON-ZEROS INPUT. IT IS NOT */ /* ALTERED. */ /* IRN(I),I=1,2,...,NZ MUST BE SET TO THE ROW INDICES OF THE */ /* NON-ZEROS ON INPUT. IT IS NOT ALTERED UNLESS EQUIVALENCED WITH IW. */ /* IRN(1) MAY BE EQUIVALENCED WITH IW(1). */ /* ICN(I),I=1,2,...,NZ MUST BE SET TO THE COLUMN INDICES OF THE */ /* NON-ZEROS ON INPUT. IT IS NOT ALTERED UNLESS EQUIVALENCED */ /* WITH IW.ICN(1) MAY BE EQUIVELENCED WITH IW(K),K.GT.NZ. */ /* PERM(I) MUST BE SET TO HOLD THE POSITION OF VARIABLE I IN THE */ /* PERMUTED ORDER. IT IS NOT ALTERED. */ /* IW NEED NOT BE SET ON INPUT. ON OUTPUT IT CONTAINS LISTS OF */ /* COLUMN NUMBERS, EACH LIST BEING HEADED BY ITS LENGTH. */ /* LW MUST BE SET TO THE LENGTH OF IW. IT MUST BE AT LEAST */ /* MAX(NZ,N+(NO. OF OFF-DIAGONAL NON-ZEROS)). IT IS NOT ALTERED. */ /* IPE NEED NOT BE SET ON INPUT. ON OUTPUT IPE(I) POINTS TO THE START OF */ /* THE ENTRY IN IW FOR ROW I, OR IS ZERO IF THERE IS NO ENTRY. */ /* IQ NEED NOT BE SET. ON OUTPUT IQ(I) CONTAINS THE NUMBER OF */ /* OFF-DIAGONAL NON-ZEROS IN ROW I, INCLUDING DUPLICATES. */ /* FLAG IS USED FOR WORKSPACE TO HOLD FLAGS TO PERMIT DUPLICATE */ /* ENTRIES TO BE IDENTIFIED QUICKLY. */ /* IWFR NEED NOT BE SET ON INPUT. ON OUTPUT IT POINTS TO THE FIRST */ /* UNUSED LOCATION IN IW. */ /* ICNTL is an INTEGER array of length 30, see MA27A/AD. */ /* INFO is an INTEGER array of length 20, see MA27A/AD. */ /* .. Scalar Arguments .. */ /* .. */ /* .. Array Arguments .. */ /* .. */ /* .. Local Scalars .. */ /* .. */ /* .. Intrinsic Functions .. */ /* .. */ /* .. Executable Statements .. */ /* INITIALIZE INFO(1), INFO(2) AND IQ */ /* Parameter adjustments */ --flag__; --iq; --ipe; --perm; --irn; --icn; --iw; --icntl; --info; /* Function Body */ info[1] = 0; info[2] = 0; i__1 = *n; for (i__ = 1; i__ <= i__1; ++i__) { iq[i__] = 0; /* L10: */ } /* COUNT THE NUMBERS OF NON-ZEROS IN THE ROWS, PRINT WARNINGS ABOUT */ /* OUT-OF-RANGE INDICES AND TRANSFER GENUINE ROW NUMBERS */ /* (NEGATED) INTO IW. */ if (*nz == 0) { goto L110; } i__1 = *nz; for (k = 1; k <= i__1; ++k) { i__ = irn[k]; j = icn[k]; iw[k] = -i__; if (i__ < j) { if (i__ >= 1 && j <= *n) { goto L80; } } else if (i__ > j) { if (j >= 1 && i__ <= *n) { goto L80; } } else { iw[k] = 0; if (i__ >= 1 && i__ <= *n) { goto L100; } } ++info[2]; info[1] = 1; iw[k] = 0; if (info[2] <= 1 && icntl[2] > 0) { io___144.ciunit = icntl[2]; s_wsfe(&io___144); do_fio(&c__1, (char *)&info[1], (ftnlen)sizeof(integer)); e_wsfe(); } if (info[2] <= 10 && icntl[2] > 0) { io___145.ciunit = icntl[2]; s_wsfe(&io___145); do_fio(&c__1, (char *)&k, (ftnlen)sizeof(integer)); do_fio(&c__1, (char *)&i__, (ftnlen)sizeof(integer)); do_fio(&c__1, (char *)&j, (ftnlen)sizeof(integer)); e_wsfe(); } goto L100; L80: if (perm[j] > perm[i__]) { goto L90; } ++iq[j]; goto L100; L90: ++iq[i__]; L100: ; } /* ACCUMULATE ROW COUNTS TO GET POINTERS TO ROW ENDS */ /* IN IPE. */ L110: *iwfr = 1; lbig = 0; i__1 = *n; for (i__ = 1; i__ <= i__1; ++i__) { l = iq[i__]; lbig = max(l,lbig); *iwfr += l; ipe[i__] = *iwfr - 1; /* L120: */ } /* PERFORM IN-PLACE SORT */ if (*nz == 0) { goto L250; } i__1 = *nz; for (k = 1; k <= i__1; ++k) { i__ = -iw[k]; if (i__ <= 0) { goto L160; } l = k; iw[k] = 0; i__2 = *nz; for (id = 1; id <= i__2; ++id) { j = icn[l]; if (perm[i__] < perm[j]) { goto L130; } l = ipe[j]; ipe[j] = l - 1; in = iw[l]; iw[l] = i__; goto L140; L130: l = ipe[i__]; ipe[i__] = l - 1; in = iw[l]; iw[l] = j; L140: i__ = -in; if (i__ <= 0) { goto L160; } /* L150: */ } L160: ; } /* MAKE ROOM IN IW FOR ROW LENGTHS AND INITIALIZE FLAG. */ k = *iwfr - 1; l = k + *n; *iwfr = l + 1; i__1 = *n; for (i__ = 1; i__ <= i__1; ++i__) { flag__[i__] = 0; j = *n + 1 - i__; len = iq[j]; if (len <= 0) { goto L180; } i__2 = len; for (jdummy = 1; jdummy <= i__2; ++jdummy) { iw[l] = iw[k]; --k; --l; /* L170: */ } L180: ipe[j] = l; --l; /* L190: */ } if (lbig >= icntl[4]) { goto L210; } /* PLACE ROW LENGTHS IN IW */ i__1 = *n; for (i__ = 1; i__ <= i__1; ++i__) { k = ipe[i__]; iw[k] = iq[i__]; if (iq[i__] == 0) { ipe[i__] = 0; } /* L200: */ } goto L250; /* REMOVE DUPLICATE ENTRIES */ L210: *iwfr = 1; i__1 = *n; for (i__ = 1; i__ <= i__1; ++i__) { k1 = ipe[i__] + 1; k2 = ipe[i__] + iq[i__]; if (k1 <= k2) { goto L220; } ipe[i__] = 0; goto L240; L220: ipe[i__] = *iwfr; ++(*iwfr); i__2 = k2; for (k = k1; k <= i__2; ++k) { j = iw[k]; if (flag__[j] == i__) { goto L230; } iw[*iwfr] = j; ++(*iwfr); flag__[j] = i__; L230: ; } k = ipe[i__]; iw[k] = *iwfr - k - 1; L240: ; } L250: return 0; } /* ma27j_ */ /* Subroutine */ int ma27k_(integer *n, integer *ipe, integer *iw, integer * lw, integer *iwfr, integer *ips, integer *ipv, integer *nv, integer * flag__, integer *ncmpa) { /* System generated locals */ integer i__1, i__2, i__3, i__4, i__5; /* Local variables */ static integer i__, j, ie, je, me, ip, jp, ln, ml, js, ms, jp1, jp2; extern /* Subroutine */ int ma27u_(integer *, integer *, integer *, integer *, integer *, integer *); static integer lwfr, minjs, kdummy; /* USING A GIVEN PIVOTAL SEQUENCE AND A REPRESENTATION OF THE MATRIX THAT */ /* INCLUDES ONLY NON-ZEROS OF THE STRICTLY UPPER-TRIANGULAR PART */ /* OF THE PERMUTED MATRIX, CONSTRUCT TREE POINTERS. */ /* N MUST BE SET TO THE MATRIX ORDER. IT IS NOT ALTERED. */ /* IPE(I) MUST BE SET TO POINT TO THE POSITION IN IW OF THE */ /* START OF ROW I OR HAVE THE VALUE ZERO IF ROW I HAS NO OFF- */ /* DIAGONAL NON-ZEROS. DURING EXECUTION IT IS USED AS FOLLOWS. */ /* IF VARIABLE I IS ELIMINATED THEN IPE(I) POINTS TO THE LIST */ /* OF VARIABLES FOR CREATED ELEMENT I. IF ELEMENT I IS */ /* ABSORBED INTO NEWLY CREATED ELEMENT J THEN IPE(I)=-J. */ /* IW MUST BE SET ON ENTRY TO HOLD LISTS OF VARIABLES BY */ /* ROWS, EACH LIST BEING HEADED BY ITS LENGTH. WHEN A VARIABLE */ /* IS ELIMINATED ITS LIST IS REPLACED BY A LIST OF VARIABLES */ /* IN THE NEW ELEMENT. */ /* LW MUST BE SET TO THE LENGTH OF IW. IT IS NOT ALTERED. */ /* IWFR MUST BE SET TO THE POSITION IN IW OF THE FIRST FREE VARIABLE. */ /* IT IS REVISED DURING EXECUTION, CONTINUING TO HAVE THIS MEANING. */ /* IPS(I) MUST BE SET TO THE POSITION OF VARIABLE I IN THE REQUIRED */ /* ORDERING. IT IS NOT ALTERED. */ /* IPV NEED NOT BE SET. IPV(K) IS SET TO HOLD THE K TH VARIABLE */ /* IN PIVOT ORDER. */ /* NV NEED NOT BE SET. IF VARIABLE J HAS NOT BEEN ELIMINATED THEN */ /* THE LAST ELEMENT WHOSE LEADING VARIABLE (VARIABLE EARLIEST */ /* IN THE PIVOT SEQUENCE) IS J IS ELEMENT NV(J). IF ELEMENT J */ /* EXISTS THEN THE LAST ELEMENT HAVING THE SAME LEADING */ /* VARIABLE IS NV(J). IN BOTH CASES NV(J)=0 IF THERE IS NO SUCH */ /* ELEMENT. IF ELEMENT J HAS BEEN MERGED INTO A LATER ELEMENT */ /* THEN NV(J) IS THE DEGREE AT THE TIME OF ELIMINATION. */ /* FLAG IS USED AS WORKSPACE FOR VARIABLE FLAGS. */ /* FLAG(JS)=ME IF JS HAS BEEN INCLUDED IN THE LIST FOR ME. */ /* NCMPA see INFO(11) in MA27A/AD. */ /* .. Scalar Arguments .. */ /* .. */ /* .. Array Arguments .. */ /* .. */ /* .. Local Scalars .. */ /* .. */ /* .. External Subroutines .. */ /* .. */ /* .. Intrinsic Functions .. */ /* .. */ /* .. Executable Statements .. */ /* INITIALIZATIONS */ /* Parameter adjustments */ --flag__; --nv; --ipv; --ips; --ipe; --iw; /* Function Body */ i__1 = *n; for (i__ = 1; i__ <= i__1; ++i__) { flag__[i__] = 0; nv[i__] = 0; j = ips[i__]; ipv[j] = i__; /* L10: */ } *ncmpa = 0; /* START OF MAIN LOOP */ i__1 = *n; for (ml = 1; ml <= i__1; ++ml) { /* ME=MS IS THE NAME OF THE VARIABLE ELIMINATED AND */ /* OF THE ELEMENT CREATED IN THE MAIN LOOP. */ ms = ipv[ml]; me = ms; flag__[ms] = me; /* MERGE ROW MS WITH ALL THE ELEMENTS HAVING MS AS LEADING VARIABLE. */ /* IP POINTS TO THE START OF THE NEW LIST. */ ip = *iwfr; /* MINJS IS SET TO THE POSITION IN THE ORDER OF THE LEADING VARIABLE */ /* IN THE NEW LIST. */ minjs = *n; ie = me; i__2 = *n; for (kdummy = 1; kdummy <= i__2; ++kdummy) { /* SEARCH VARIABLE LIST OF ELEMENT IE. */ /* JP POINTS TO THE CURRENT POSITION IN THE LIST BEING SEARCHED. */ jp = ipe[ie]; /* LN IS THE LENGTH OF THE LIST BEING SEARCHED. */ ln = 0; if (jp <= 0) { goto L60; } ln = iw[jp]; /* SEARCH FOR DIFFERENT VARIABLES AND ADD THEM TO LIST, */ /* COMPRESSING WHEN NECESSARY */ i__3 = ln; for (jp1 = 1; jp1 <= i__3; ++jp1) { ++jp; /* PLACE NEXT VARIABLE IN JS. */ js = iw[jp]; /* JUMP IF VARIABLE HAS ALREADY BEEN INCLUDED. */ if (flag__[js] == me) { goto L50; } flag__[js] = me; if (*iwfr < *lw) { goto L40; } /* PREPARE FOR COMPRESSING IW BY ADJUSTING POINTER TO AND LENGTH OF */ /* THE LIST FOR IE TO REFER TO THE REMAINING ENTRIES. */ ipe[ie] = jp; iw[jp] = ln - jp1; /* COMPRESS IW. */ i__4 = ip - 1; ma27u_(n, &ipe[1], &iw[1], &i__4, &lwfr, ncmpa); /* COPY NEW LIST FORWARD */ jp2 = *iwfr - 1; *iwfr = lwfr; if (ip > jp2) { goto L30; } i__4 = jp2; for (jp = ip; jp <= i__4; ++jp) { iw[*iwfr] = iw[jp]; ++(*iwfr); /* L20: */ } L30: ip = lwfr; jp = ipe[ie]; /* ADD VARIABLE JS TO NEW LIST. */ L40: iw[*iwfr] = js; /* Computing MIN */ i__4 = minjs, i__5 = ips[js]; minjs = min(i__4,i__5); ++(*iwfr); L50: ; } /* RECORD ABSORPTION OF ELEMENT IE INTO NEW ELEMENT. */ L60: ipe[ie] = -me; /* PICK UP NEXT ELEMENT WITH LEADING VARIABLE MS. */ je = nv[ie]; /* STORE DEGREE OF IE. */ nv[ie] = ln + 1; ie = je; /* LEAVE LOOP IF THERE ARE NO MORE ELEMENTS. */ if (ie == 0) { goto L80; } /* L70: */ } L80: if (*iwfr > ip) { goto L90; } /* DEAL WITH NULL NEW ELEMENT. */ ipe[me] = 0; nv[me] = 1; goto L100; /* LINK NEW ELEMENT WITH OTHERS HAVING SAME LEADING VARIABLE. */ L90: minjs = ipv[minjs]; nv[me] = nv[minjs]; nv[minjs] = me; /* MOVE FIRST ENTRY IN NEW LIST TO END TO ALLOW ROOM FOR LENGTH AT */ /* FRONT. SET POINTER TO FRONT. */ iw[*iwfr] = iw[ip]; iw[ip] = *iwfr - ip; ipe[me] = ip; ++(*iwfr); L100: ; } return 0; } /* ma27k_ */ /* Subroutine */ int ma27l_(integer *n, integer *ipe, integer *nv, integer * ips, integer *ne, integer *na, integer *nd, integer *nsteps, integer * nemin) { /* System generated locals */ integer i__1, i__2; /* Local variables */ static integer i__, k, l, ib, if__, il, is, nr, ison; /* TREE SEARCH */ /* GIVEN SON TO FATHER TREE POINTERS, PERFORM DEPTH-FIRST */ /* SEARCH TO FIND PIVOT ORDER AND NUMBER OF ELIMINATIONS */ /* AND ASSEMBLIES AT EACH STAGE. */ /* N MUST BE SET TO THE MATRIX ORDER. IT IS NOT ALTERED. */ /* IPE(I) MUST BE SET EQUAL TO -(FATHER OF NODE I) OR ZERO IF */ /* NODE IS A ROOT. IT IS ALTERED TO POINT TO ITS NEXT */ /* YOUNGER BROTHER IF IT HAS ONE, BUT OTHERWISE IS NOT */ /* CHANGED. */ /* NV(I) MUST BE SET TO ZERO IF NO VARIABLES ARE ELIMINATED AT NODE */ /* I AND TO THE DEGREE OTHERWISE. ONLY LEAF NODES CAN HAVE */ /* ZERO VALUES OF NV(I). NV IS NOT ALTERED. */ /* IPS(I) NEED NOT BE SET. IT IS USED TEMPORARILY TO HOLD */ /* -(ELDEST SON OF NODE I) IF IT HAS ONE AND 0 OTHERWISE. IT IS */ /* EVENTUALLY SET TO HOLD THE POSITION OF NODE I IN THE ORDER. */ /* NE(IS) NEED NOT BE SET. IT IS SET TO THE NUMBER OF VARIABLES */ /* ELIMINATED AT STAGE IS OF THE ELIMINATION. */ /* NA(IS) NEED NOT BE SET. IT IS SET TO THE NUMBER OF ELEMENTS */ /* ASSEMBLED AT STAGE IS OF THE ELIMINATION. */ /* ND(IS) NEED NOT BE SET. IT IS SET TO THE DEGREE AT STAGE IS OF */ /* THE ELIMINATION. */ /* NSTEPS NEED NOT BE SET. IT IS SET TO THE NUMBER OF ELIMINATION */ /* STEPS. */ /* NEMIN see ICNTL(5) in MA27A/AD. */ /* .. Scalar Arguments .. */ /* .. */ /* .. Array Arguments .. */ /* .. */ /* .. Local Scalars .. */ /* .. */ /* .. Executable Statements .. */ /* INITIALIZE IPS AND NE. */ /* Parameter adjustments */ --nd; --na; --ne; --ips; --nv; --ipe; /* Function Body */ i__1 = *n; for (i__ = 1; i__ <= i__1; ++i__) { ips[i__] = 0; ne[i__] = 0; /* L10: */ } /* SET IPS(I) TO -(ELDEST SON OF NODE I) AND IPE(I) TO NEXT YOUNGER */ /* BROTHER OF NODE I IF IT HAS ONE. */ /* FIRST PASS IS FOR NODES WITHOUT ELIMINATIONS. */ i__1 = *n; for (i__ = 1; i__ <= i__1; ++i__) { if (nv[i__] > 0) { goto L20; } if__ = -ipe[i__]; is = -ips[if__]; if (is > 0) { ipe[i__] = is; } ips[if__] = -i__; L20: ; } /* NR IS DECREMENTED FOR EACH ROOT NODE. THESE ARE STORED IN */ /* NE(I),I=NR,N. */ nr = *n + 1; /* SECOND PASS TO ADD NODES WITH ELIMINATIONS. */ i__1 = *n; for (i__ = 1; i__ <= i__1; ++i__) { if (nv[i__] <= 0) { goto L50; } /* NODE IF IS THE FATHER OF NODE I. */ if__ = -ipe[i__]; if (if__ == 0) { goto L40; } is = -ips[if__]; /* JUMP IF NODE IF HAS NO SONS YET. */ if (is <= 0) { goto L30; } /* SET POINTER TO NEXT BROTHER */ ipe[i__] = is; /* NODE I IS ELDEST SON OF NODE IF. */ L30: ips[if__] = -i__; goto L50; /* WE HAVE A ROOT */ L40: --nr; ne[nr] = i__; L50: ; } /* DEPTH-FIRST SEARCH. */ /* IL HOLDS THE CURRENT TREE LEVEL. ROOTS ARE AT LEVEL N, THEIR SONS */ /* ARE AT LEVEL N-1, ETC. */ /* IS HOLDS THE CURRENT ELIMINATION STAGE. WE ACCUMULATE THE NUMBER */ /* OF ELIMINATIONS AT STAGE IS DIRECTLY IN NE(IS). THE NUMBER OF */ /* ASSEMBLIES IS ACCUMULATED TEMPORARILY IN NA(IL), FOR TREE */ /* LEVEL IL, AND IS TRANSFERED TO NA(IS) WHEN WE REACH THE */ /* APPROPRIATE STAGE IS. */ is = 1; /* I IS THE CURRENT NODE. */ i__ = 0; i__1 = *n; for (k = 1; k <= i__1; ++k) { if (i__ > 0) { goto L60; } /* PICK UP NEXT ROOT. */ i__ = ne[nr]; ne[nr] = 0; ++nr; il = *n; na[*n] = 0; /* GO TO SON FOR AS LONG AS POSSIBLE, CLEARING FATHER-SON POINTERS */ /* IN IPS AS EACH IS USED AND SETTING NA(IL)=0 FOR ALL LEVELS */ /* REACHED. */ L60: i__2 = *n; for (l = 1; l <= i__2; ++l) { if (ips[i__] >= 0) { goto L80; } ison = -ips[i__]; ips[i__] = 0; i__ = ison; --il; na[il] = 0; /* L70: */ } /* RECORD POSITION OF NODE I IN THE ORDER. */ L80: ips[i__] = k; ++ne[is]; /* JUMP IF NODE HAS NO ELIMINATIONS. */ if (nv[i__] <= 0) { goto L120; } if (il < *n) { ++na[il + 1]; } na[is] = na[il]; nd[is] = nv[i__]; /* CHECK FOR STATIC CONDENSATION */ if (na[is] != 1) { goto L90; } if (nd[is - 1] - ne[is - 1] == nd[is]) { goto L100; } /* CHECK FOR SMALL NUMBERS OF ELIMINATIONS IN BOTH LAST TWO STEPS. */ L90: if (ne[is] >= *nemin) { goto L110; } if (na[is] == 0) { goto L110; } if (ne[is - 1] >= *nemin) { goto L110; } /* COMBINE THE LAST TWO STEPS */ L100: na[is - 1] = na[is - 1] + na[is] - 1; nd[is - 1] = nd[is] + ne[is - 1]; ne[is - 1] = ne[is] + ne[is - 1]; ne[is] = 0; goto L120; L110: ++is; L120: ib = ipe[i__]; if (ib >= 0) { /* NODE I HAS A BROTHER OR IS A ROOT */ if (ib > 0) { na[il] = 0; } i__ = ib; } else { /* GO TO FATHER OF NODE I */ i__ = -ib; ++il; } /* L160: */ } *nsteps = is - 1; return 0; } /* ma27l_ */ /* Subroutine */ int ma27m_(integer *n, integer *nz, integer *irn, integer * icn, integer *perm, integer *na, integer *ne, integer *nd, integer * nsteps, integer *lstki, integer *lstkr, integer *iw, integer *info, real *ops) { /* System generated locals */ integer i__1, i__2, i__3; /* Local variables */ static integer i__, k, nz1, nz2, nfr, iold, jold, iorg, jorg, inew, itop, lstk, nstk, irow; static real delim; static integer nelim, itree, istki, nassr, istkr, nirnec, nrlnec, niradu, nrladu, numorg, nirtot, nrltot; /* STORAGE AND OPERATION COUNT EVALUATION. */ /* EVALUATE NUMBER OF OPERATIONS AND SPACE REQUIRED BY FACTORIZATION */ /* USING MA27B/BD. THE VALUES GIVEN ARE EXACT ONLY IF NO NUMERICAL */ /* PIVOTING IS PERFORMED AND THEN ONLY IF IRN(1) WAS NOT */ /* EQUIVALENCED TO IW(1) BY THE USER BEFORE CALLING MA27A/AD. IF */ /* THE EQUIVALENCE HAS BEEN MADE ONLY AN UPPER BOUND FOR NIRNEC */ /* AND NRLNEC CAN BE CALCULATED ALTHOUGH THE OTHER COUNTS WILL */ /* STILL BE EXACT. */ /* N MUST BE SET TO THE MATRIX ORDER. IT IS NOT ALTERED. */ /* NZ MUST BE SET TO THE NUMBER OF NON-ZEROS INPUT. IT IS NOT ALTERED. */ /* IRN,ICN. UNLESS IRN(1) HAS BEEN EQUIVALENCED TO IW(1) */ /* IRN,ICN MUST BE SET TO THE ROW AND COLUMN INDICES OF THE */ /* NON-ZEROS ON INPUT. THEY ARE NOT ALTERED BY MA27M/MD. */ /* PERM MUST BE SET TO THE POSITION IN THE PIVOT ORDER OF EACH ROW. */ /* IT IS NOT ALTERED. */ /* NA,NE,ND MUST BE SET TO HOLD, FOR EACH TREE NODE, THE NUMBER OF STACK */ /* ELEMENTS ASSEMBLED, THE NUMBER OF ELIMINATIONS AND THE SIZE OF */ /* THE ASSEMBLED FRONT MATRIX RESPECTIVELY. THEY ARE NOT ALTERED. */ /* NSTEPS MUST BE SET TO HOLD THE NUMBER OF TREE NODES. IT IS NOT */ /* ALTERED. */ /* LSTKI IS USED AS A WORK ARRAY BY MA27M/MD. */ /* LSTKR. IF IRN(1) IS EQUIVALENCED TO IW(1) THEN LSTKR(I) */ /* MUST HOLD THE TOTAL NUMBER OF OFF-DIAGONAL ENTRIES (INCLUDING */ /* DUPLICATES) IN ROW I (I=1,..,N) OF THE ORIGINAL MATRIX. IT */ /* IS USED AS WORKSPACE BY MA27M/MD. */ /* IW IS A WORKSPACE ARRAY USED BY OTHER SUBROUTINES AND PASSED TO THIS */ /* SUBROUTINE ONLY SO THAT A TEST FOR EQUIVALENCE WITH IRN CAN BE */ /* MADE. */ /* COUNTS FOR OPERATIONS AND STORAGE ARE ACCUMULATED IN VARIABLES */ /* OPS,NRLTOT,NIRTOT,NRLNEC,NIRNEC,NRLADU,NRLNEC,NIRNEC. */ /* OPS NUMBER OF MULTIPLICATIONS AND ADDITIONS DURING FACTORIZATION. */ /* NRLADU,NIRADU REAL AND INTEGER STORAGE RESPECTIVELY FOR THE */ /* MATRIX FACTORS. */ /* NRLTOT,NIRTOT REAL AND INTEGER STRORAGE RESPECTIVELY REQUIRED */ /* FOR THE FACTORIZATION IF NO COMPRESSES ARE ALLOWED. */ /* NRLNEC,NIRNEC REAL AND INTEGER STORAGE RESPECTIVELY REQUIRED FOR */ /* THE FACTORIZATION IF COMPRESSES ARE ALLOWED. */ /* INFO is an INTEGER array of length 20, see MA27A/AD. */ /* OPS ACCUMULATES THE NO. OF MULTIPLY/ADD PAIRS NEEDED TO CREATE THE */ /* TRIANGULAR FACTORIZATION, IN THE DEFINITE CASE. */ /* .. Scalar Arguments .. */ /* .. */ /* .. Array Arguments .. */ /* .. */ /* .. Local Scalars .. */ /* .. */ /* .. Intrinsic Functions .. */ /* .. */ /* .. Executable Statements .. */ /* Parameter adjustments */ --lstkr; --lstki; --perm; --irn; --icn; --nd; --ne; --na; --iw; --info; /* Function Body */ if (*nz == 0) { goto L20; } /* JUMP IF IW AND IRN HAVE NOT BEEN EQUIVALENCED. */ if (irn[1] != iw[1]) { goto L20; } /* RESET IRN(1). */ irn[1] = iw[1] - 1; /* THE TOTAL NUMBER OF OFF-DIAGONAL ENTRIES IS ACCUMULATED IN NZ2. */ /* LSTKI IS SET TO HOLD THE TOTAL NUMBER OF ENTRIES (INCUDING */ /* THE DIAGONAL) IN EACH ROW IN PERMUTED ORDER. */ nz2 = 0; i__1 = *n; for (iold = 1; iold <= i__1; ++iold) { inew = perm[iold]; lstki[inew] = lstkr[iold] + 1; nz2 += lstkr[iold]; /* L10: */ } /* NZ1 IS THE NUMBER OF ENTRIES IN ONE TRIANGLE INCLUDING THE DIAGONAL. */ /* NZ2 IS THE TOTAL NUMBER OF ENTRIES INCLUDING THE DIAGONAL. */ nz1 = nz2 / 2 + *n; nz2 += *n; goto L60; /* COUNT (IN LSTKI) NON-ZEROS IN ORIGINAL MATRIX BY PERMUTED ROW (UPPER */ /* TRIANGLE ONLY). INITIALIZE COUNTS. */ L20: i__1 = *n; for (i__ = 1; i__ <= i__1; ++i__) { lstki[i__] = 1; /* L30: */ } /* ACCUMULATE NUMBER OF NON-ZEROS WITH INDICES IN RANGE IN NZ1 */ /* DUPLICATES ON THE DIAGONAL ARE IGNORED BUT NZ1 INCLUDES ANY */ /* DIAGONALS NOT PRESENT ON INPUT. */ /* ACCUMULATE ROW COUNTS IN LSTKI. */ nz1 = *n; if (*nz == 0) { goto L50; } i__1 = *nz; for (i__ = 1; i__ <= i__1; ++i__) { iold = irn[i__]; jold = icn[i__]; /* JUMP IF INDEX IS OUT OF RANGE. */ if (iold < 1 || iold > *n) { goto L40; } if (jold < 1 || jold > *n) { goto L40; } if (iold == jold) { goto L40; } ++nz1; /* Computing MIN */ i__2 = perm[iold], i__3 = perm[jold]; irow = min(i__2,i__3); ++lstki[irow]; L40: ; } L50: nz2 = nz1; /* ISTKR,ISTKI CURRENT NUMBER OF STACK ENTRIES IN */ /* REAL AND INTEGER STORAGE RESPECTIVELY. */ /* OPS,NRLADU,NIRADU,NIRTOT,NRLTOT,NIRNEC,NRLNEC,NZ2 ARE DEFINED ABOVE. */ /* NZ2 CURRENT NUMBER OF ORIGINAL MATRIX ENTRIES NOT YET PROCESSED. */ /* NUMORG CURRENT TOTAL NUMBER OF ROWS ELIMINATED. */ /* ITOP CURRENT NUMBER OF ELEMENTS ON THE STACK. */ L60: istki = 0; istkr = 0; *ops = 0.f; nrladu = 0; /* ONE LOCATION IS NEEDED TO RECORD THE NUMBER OF BLOCKS */ /* ACTUALLY USED. */ niradu = 1; nirtot = nz1; nrltot = nz1; nirnec = nz2; nrlnec = nz2; numorg = 0; itop = 0; /* EACH PASS THROUGH THIS LOOP PROCESSES A NODE OF THE TREE. */ i__1 = *nsteps; for (itree = 1; itree <= i__1; ++itree) { nelim = ne[itree]; delim = (real) nelim; nfr = nd[itree]; nstk = na[itree]; /* ADJUST STORAGE COUNTS ON ASSEMBLY OF CURRENT FRONTAL MATRIX. */ nassr = nfr * (nfr + 1) / 2; if (nstk != 0) { nassr = nassr - lstkr[itop] + 1; } /* Computing MAX */ i__2 = nrltot, i__3 = nrladu + nassr + istkr + nz1; nrltot = max(i__2,i__3); /* Computing MAX */ i__2 = nirtot, i__3 = niradu + nfr + 2 + istki + nz1; nirtot = max(i__2,i__3); /* Computing MAX */ i__2 = nrlnec, i__3 = nrladu + nassr + istkr + nz2; nrlnec = max(i__2,i__3); /* Computing MAX */ i__2 = nirnec, i__3 = niradu + nfr + 2 + istki + nz2; nirnec = max(i__2,i__3); /* DECREASE NZ2 BY THE NUMBER OF ENTRIES IN ROWS BEING ELIMINATED AT */ /* THIS STAGE. */ i__2 = nelim; for (iorg = 1; iorg <= i__2; ++iorg) { jorg = numorg + iorg; nz2 -= lstki[jorg]; /* L70: */ } numorg += nelim; /* JUMP IF THERE ARE NO STACK ASSEMBLIES AT THIS NODE. */ if (nstk <= 0) { goto L90; } /* REMOVE ELEMENTS FROM THE STACK. THERE ARE ITOP ELEMENTS ON THE */ /* STACK WITH THE APPROPRIATE ENTRIES IN LSTKR,LSTKI GIVING */ /* THE REAL AND INTEGER STORAGE RESPECTIVELY FOR EACH STACK */ /* ELEMENT. */ i__2 = nstk; for (k = 1; k <= i__2; ++k) { lstk = lstkr[itop]; istkr -= lstk; lstk = lstki[itop]; istki -= lstk; --itop; /* L80: */ } /* ACCUMULATE NON-ZEROS IN FACTORS AND NUMBER OF OPERATIONS. */ L90: nrladu += nelim * ((nfr << 1) - nelim + 1) / 2; niradu = niradu + 2 + nfr; if (nelim == 1) { --niradu; } *ops += (nfr * delim * (nfr + 1) - ((nfr << 1) + 1) * delim * (delim + 1) / 2 + delim * (delim + 1) * (delim * 2 + 1) / 6) / 2; if (itree == *nsteps) { goto L100; } /* JUMP IF ALL OF FRONTAL MATRIX HAS BEEN ELIMINATED. */ if (nfr == nelim) { goto L100; } /* STACK REMAINDER OF ELEMENT. */ ++itop; lstkr[itop] = (nfr - nelim) * (nfr - nelim + 1) / 2; lstki[itop] = nfr - nelim + 1; istki += lstki[itop]; istkr += lstkr[itop]; /* WE DO NOT NEED TO ADJUST THE COUNTS FOR THE REAL STORAGE BECAUSE */ /* THE REMAINDER OF THE FRONTAL MATRIX IS SIMPLY MOVED IN THE */ /* STORAGE FROM FACTORS TO STACK AND NO EXTRA STORAGE IS REQUIRED. */ /* Computing MAX */ i__2 = nirtot, i__3 = niradu + istki + nz1; nirtot = max(i__2,i__3); /* Computing MAX */ i__2 = nirnec, i__3 = niradu + istki + nz2; nirnec = max(i__2,i__3); L100: ; } /* ADJUST THE STORAGE COUNTS TO ALLOW FOR THE USE OF THE REAL AND */ /* INTEGER STORAGE FOR PURPOSES OTHER THAN PURELY THE */ /* FACTORIZATION ITSELF. */ /* THE SECOND TWO TERMS ARE THE MINUMUM AMOUNT REQUIRED BY MA27N/ND. */ /* Computing MAX */ i__1 = nrlnec, i__2 = *n + max(*nz,nz1); nrlnec = max(i__1,i__2); /* Computing MAX */ i__1 = nrltot, i__2 = *n + max(*nz,nz1); nrltot = max(i__1,i__2); nrlnec = min(nrlnec,nrltot); nirnec = max(*nz,nirnec); nirtot = max(*nz,nirtot); nirnec = min(nirnec,nirtot); info[3] = nrltot; info[4] = nirtot; info[5] = nrlnec; info[6] = nirnec; info[7] = nrladu; info[8] = niradu; return 0; } /* ma27m_ */ /* Subroutine */ int ma27n_(integer *n, integer *nz, integer *nz1, real *a, integer *la, integer *irn, integer *icn, integer *iw, integer *liw, integer *perm, integer *iw2, integer *icntl, integer *info) { /* Format strings */ static char fmt_40[] = "(\002 *** WARNING MESSAGE FROM SUBROUTINE MA27" "B\002,\002 *** INFO(1) =\002,i2)"; static char fmt_50[] = "(i6,\002TH NON-ZERO (IN ROW\002,i6,\002 AND COLU" "MN\002,i6,\002) IGNORED\002)"; /* System generated locals */ integer i__1, i__2; /* Builtin functions */ integer s_wsfe(cilist *), do_fio(integer *, char *, ftnlen), e_wsfe(void); /* Local variables */ static integer i__, k, j1, j2, ia, ii, jj, ich, iiw, iold, jold, inew; static real anow; static integer jnew, ipos, jpos; static real anext; /* Fortran I/O blocks */ static cilist io___212 = { 0, 0, 0, fmt_40, 0 }; static cilist io___213 = { 0, 0, 0, fmt_50, 0 }; /* SORT PRIOR TO FACTORIZATION USING MA27O/OD. */ /* THIS SUBROUTINE REORDERS THE USER'S INPUT SO THAT THE UPPER TRIANGLE */ /* OF THE PERMUTED MATRIX, INCLUDING THE DIAGONAL, IS */ /* HELD ORDERED BY ROWS AT THE END OF THE STORAGE FOR A AND IW. */ /* IT IGNORES ENTRIES WITH ONE OR BOTH INDICES OUT OF RANGE AND */ /* ACCUMULATES DIAGONAL ENTRIES. */ /* IT ADDS EXPLICIT ZEROS ON THE DIAGONAL WHERE NECESSARY. */ /* N - MUST BE SET TO THE ORDER OF THE MATRIX. */ /* IT IS NOT ALTERED BY MA27N/ND. */ /* NZ - ON ENTRY NZ MUST BE SET TO THE NUMBER */ /* OF NON-ZEROS INPUT. NOT ALTERED BY MA27N/ND. */ /* NZ1 - ON EXIT NZ1 WILL BE EQUAL TO THE NUMBER OF ENTRIES IN THE */ /* SORTED MATRIX. */ /* A - ON ENTRY A(I) MUST */ /* HOLD THE VALUE OF THE ORIGINAL MATRIX ELEMENT IN POSITION */ /* (IRN(I),ICN(I)),I=1,NZ. ON EXIT A(LA-NZ1+I),I=1,NZ1 HOLDS */ /* THE UPPER TRIANGLE OF THE PERMUTED MATRIX BY ROWS WITH */ /* THE DIAGONAL ENTRY FIRST ALTHOUGH THERE IS NO FURTHER */ /* ORDERING WITHIN THE ROWS THEMSELVES. */ /* LA - LENGTH OF ARRAY A. MUST BE AT LEAST N+MAX(NZ,NZ1). */ /* IT IS NOT ALTERED BY MA27N/ND. */ /* IRN - IRN(I) MUST BE SET TO */ /* HOLD THE ROW INDEX OF ENTRY A(I),I=1,NZ. IRN WILL BE */ /* UNALTERED BY MA27N/ND, UNLESS IT IS EQUIVALENCED WITH IW. */ /* ICN - ICN(I) MUST BE SET TO */ /* HOLD THE COLUMN INDEX OF ENTRY A(I),I=1,NZ. ICN WILL BE */ /* UNALTERED BY MA27N/ND, UNLESS IT IS EQUIVALENCED WITH IW. */ /* IW - USED AS WORKSPACE AND ON */ /* EXIT, ENTRIES IW(LIW-NZ1+I),I=1,NZ1 HOLD THE COLUMN INDICES */ /* (THE ORIGINAL UNPERMUTED INDICES) OF THE CORRESPONDING ENTRY */ /* OF A WITH THE FIRST ENTRY FOR EACH ROW FLAGGED NEGATIVE. */ /* IRN(1) MAY BE EQUIVALENCED WITH IW(1) AND ICN(1) MAY BE */ /* EQUIVALENCED WITH IW(K) WHERE K.GT.NZ. */ /* LIW - LENGTH OF ARRAY IW. MUST BE AT LEAST AS */ /* GREAT AS THE MAXIMUM OF NZ AND NZ1. */ /* NOT ALTERED BY MA27N/ND. */ /* PERM - PERM(I) HOLDS THE */ /* POSITION IN THE TENTATIVE PIVOT ORDER OF ROW I IN THE */ /* ORIGINAL MATRIX,I=1,N. NOT ALTERED BY MA27N/ND. */ /* IW2 - USED AS WORKSPACE. */ /* SEE COMMENTS IN CODE IMMEDIATELY PRIOR TO */ /* EACH USE. */ /* ICNTL is an INTEGER array of length 30, see MA27A/AD. */ /* INFO is an INTEGER array of length 20, see MA27A/AD. */ /* INFO(1) - ON EXIT FROM MA27N/ND, A ZERO VALUE OF */ /* INFO(1) INDICATES THAT NO ERROR HAS BEEN DETECTED. */ /* POSSIBLE NON-ZERO VALUES ARE .. */ /* +1 WARNING. INDICES OUT OF RANGE. THESE ARE IGNORED, */ /* THEIR NUMBER IS RECORDED IN INFO(2) OF MA27E/ED AND */ /* MESSAGES IDENTIFYING THE FIRST TEN ARE OUTPUT ON UNIT */ /* ICNTL(2). */ /* -3 INTEGER ARRAY IW IS TOO SMALL. */ /* -4 REAL ARRAY A IS TOO SMALL. */ /* .. Parameters .. */ /* .. */ /* .. Scalar Arguments .. */ /* .. */ /* .. Array Arguments .. */ /* .. */ /* .. Local Scalars .. */ /* .. */ /* .. Intrinsic Functions .. */ /* .. */ /* .. Executable Statements .. */ /* Parameter adjustments */ --iw2; --perm; --a; --irn; --icn; --iw; --icntl; --info; /* Function Body */ info[1] = 0; /* INITIALIZE WORK ARRAY (IW2) IN PREPARATION FOR */ /* COUNTING NUMBERS OF NON-ZEROS IN THE ROWS AND INITIALIZE */ /* LAST N ENTRIES IN A WHICH WILL HOLD THE DIAGONAL ENTRIES */ ia = *la; i__1 = *n; for (iold = 1; iold <= i__1; ++iold) { iw2[iold] = 1; a[ia] = 0.f; --ia; /* L10: */ } /* SCAN INPUT COPYING ROW INDICES FROM IRN TO THE FIRST NZ POSITIONS */ /* IN IW. THE NEGATIVE OF THE INDEX IS HELD TO FLAG ENTRIES FOR */ /* THE IN-PLACE SORT. ENTRIES IN IW CORRESPONDING TO DIAGONALS AND */ /* ENTRIES WITH OUT-OF-RANGE INDICES ARE SET TO ZERO. */ /* FOR DIAGONAL ENTRIES, REALS ARE ACCUMULATED IN THE LAST N */ /* LOCATIONS OF A. */ /* THE NUMBER OF ENTRIES IN EACH ROW OF THE PERMUTED MATRIX IS */ /* ACCUMULATED IN IW2. */ /* INDICES OUT OF RANGE ARE IGNORED AFTER BEING COUNTED AND */ /* AFTER APPROPRIATE MESSAGES HAVE BEEN PRINTED. */ info[2] = 0; /* NZ1 IS THE NUMBER OF NON-ZEROS HELD AFTER INDICES OUT OF RANGE HAVE */ /* BEEN IGNORED AND DIAGONAL ENTRIES ACCUMULATED. */ *nz1 = *n; if (*nz == 0) { goto L80; } i__1 = *nz; for (k = 1; k <= i__1; ++k) { iold = irn[k]; if (iold > *n || iold <= 0) { goto L30; } jold = icn[k]; if (jold > *n || jold <= 0) { goto L30; } inew = perm[iold]; jnew = perm[jold]; if (inew != jnew) { goto L20; } ia = *la - *n + iold; a[ia] += a[k]; goto L60; L20: inew = min(inew,jnew); /* INCREMENT NUMBER OF ENTRIES IN ROW INEW. */ ++iw2[inew]; iw[k] = -iold; ++(*nz1); goto L70; /* ENTRY OUT OF RANGE. IT WILL BE IGNORED AND A FLAG SET. */ L30: info[1] = 1; ++info[2]; if (info[2] <= 1 && icntl[2] > 0) { io___212.ciunit = icntl[2]; s_wsfe(&io___212); do_fio(&c__1, (char *)&info[1], (ftnlen)sizeof(integer)); e_wsfe(); } if (info[2] <= 10 && icntl[2] > 0) { io___213.ciunit = icntl[2]; s_wsfe(&io___213); do_fio(&c__1, (char *)&k, (ftnlen)sizeof(integer)); do_fio(&c__1, (char *)&irn[k], (ftnlen)sizeof(integer)); do_fio(&c__1, (char *)&icn[k], (ftnlen)sizeof(integer)); e_wsfe(); } L60: iw[k] = 0; L70: ; } /* CALCULATE POINTERS (IN IW2) TO THE POSITION IMMEDIATELY AFTER THE END */ /* OF EACH ROW. */ L80: if (*nz < *nz1 && *nz1 != *n) { goto L100; } /* ROOM IS INCLUDED FOR THE DIAGONALS. */ k = 1; i__1 = *n; for (i__ = 1; i__ <= i__1; ++i__) { k += iw2[i__]; iw2[i__] = k; /* L90: */ } goto L120; /* ROOM IS NOT INCLUDED FOR THE DIAGONALS. */ L100: k = 1; i__1 = *n; for (i__ = 1; i__ <= i__1; ++i__) { k = k + iw2[i__] - 1; iw2[i__] = k; /* L110: */ } /* FAIL IF INSUFFICIENT SPACE IN ARRAYS A OR IW. */ L120: if (*nz1 > *liw) { goto L210; } if (*nz1 + *n > *la) { goto L220; } /* NOW RUN THROUGH NON-ZEROS IN ORDER PLACING THEM IN THEIR NEW */ /* POSITION AND DECREMENTING APPROPRIATE IW2 ENTRY. IF WE ARE */ /* ABOUT TO OVERWRITE AN ENTRY NOT YET MOVED, WE MUST DEAL WITH */ /* THIS AT THIS TIME. */ if (*nz1 == *n) { goto L180; } i__1 = *nz; for (k = 1; k <= i__1; ++k) { iold = -iw[k]; if (iold <= 0) { goto L140; } jold = icn[k]; anow = a[k]; iw[k] = 0; i__2 = *nz; for (ich = 1; ich <= i__2; ++ich) { inew = perm[iold]; jnew = perm[jold]; inew = min(inew,jnew); if (inew == perm[jold]) { jold = iold; } jpos = iw2[inew] - 1; iold = -iw[jpos]; anext = a[jpos]; a[jpos] = anow; iw[jpos] = jold; iw2[inew] = jpos; if (iold == 0) { goto L140; } anow = anext; jold = icn[jpos]; /* L130: */ } L140: ; } if (*nz >= *nz1) { goto L180; } /* MOVE UP ENTRIES TO ALLOW FOR DIAGONALS. */ ipos = *nz1; jpos = *nz1 - *n; i__1 = *n; for (ii = 1; ii <= i__1; ++ii) { i__ = *n - ii + 1; j1 = iw2[i__]; j2 = jpos; if (j1 > jpos) { goto L160; } i__2 = j2; for (jj = j1; jj <= i__2; ++jj) { iw[ipos] = iw[jpos]; a[ipos] = a[jpos]; --ipos; --jpos; /* L150: */ } L160: iw2[i__] = ipos + 1; --ipos; /* L170: */ } /* RUN THROUGH ROWS INSERTING DIAGONAL ENTRIES AND FLAGGING BEGINNING */ /* OF EACH ROW BY NEGATING FIRST COLUMN INDEX. */ L180: i__1 = *n; for (iold = 1; iold <= i__1; ++iold) { inew = perm[iold]; jpos = iw2[inew] - 1; ia = *la - *n + iold; a[jpos] = a[ia]; iw[jpos] = -iold; /* L190: */ } /* MOVE SORTED MATRIX TO THE END OF THE ARRAYS. */ ipos = *nz1; ia = *la; iiw = *liw; i__1 = *nz1; for (i__ = 1; i__ <= i__1; ++i__) { a[ia] = a[ipos]; iw[iiw] = iw[ipos]; --ipos; --ia; --iiw; /* L200: */ } goto L230; /* **** ERROR RETURN **** */ L210: info[1] = -3; info[2] = *nz1; goto L230; L220: info[1] = -4; info[2] = *nz1 + *n; L230: return 0; } /* ma27n_ */ /* Subroutine */ int ma27o_(integer *n, integer *nz, real *a, integer *la, integer *iw, integer *liw, integer *perm, integer *nstk, integer * nsteps, integer *maxfrt, integer *nelim, integer *iw2, integer *icntl, real *cntl, integer *info) { /* Format strings */ static char fmt_310[] = "(\002 *** WARNING MESSAGE FROM SUBROUTINE MA27" "B\002,\002 *** INFO(1) =\002,i2,/,\002 PIVOT\002,i6,\002 HAS DI" "FFERENT SIGN FROM THE PREVIOUS ONE\002)"; /* System generated locals */ integer i__1, i__2, i__3, i__4, i__5; real r__1, r__2, r__3, r__4; /* Builtin functions */ integer s_wsfe(cilist *), do_fio(integer *, char *, ftnlen), e_wsfe(void); /* Local variables */ static integer i__, j, k, j1, j2, k1, k2, jj, kk, lt; static real uu; static integer jj1, jjj, ifr, ibeg, iend, neig; static real amax; static integer iell, jcol, nblk; extern /* Subroutine */ int ma27p_(real *, integer *, integer *, integer * , integer *, integer *, integer *, integer *); static integer iass, iorg, apos, astk, jmax, jnew; static real rmax; static integer ipiv; static real tmax; static integer ipos, istk, jpiv, jpos, kmax, irow, krow, nass, npiv, ntwo; static real swop; static integer apos1, apos2, apos3, astk2, istk2, laell, iexch, liell, newel, jlast, azero, lnass; static real amult; static integer ipmnp, jnext, lnpiv, iswop, iwpos, lapos2; static real amult1, amult2; static integer npivp1, pospv1, pospv2, ncmpbi, posfac, ncmpbr, nirbdu, nrlbdu, ioldps; static real detpiv, thresh; static integer ainput, jfirst, idummy, jdummy, jmxmip, kdummy, iinput, isnpiv, nfront, numass, numorg, numstk, pivsiz, ltopst, ntotpv; /* Fortran I/O blocks */ static cilist io___280 = { 0, 0, 0, fmt_310, 0 }; /* FACTORIZATION SUBROUTINE */ /* THIS SUBROUTINE OPERATES ON THE INPUT MATRIX ORDERED BY MA27N/ND AND */ /* PRODUCES THE FACTORS OF U AND D ('A'=UTRANSPOSE*D*U) FOR USE IN */ /* SUBSEQUENT SOLUTIONS. GAUSSIAN ELIMINATION IS USED WITH PIVOTS */ /* CHOSEN FROM THE DIAGONAL. TO ENSURE STABILITY, BLOCK PIVOTS OF */ /* ORDER 2 WILL BE USED IF THE DIAGONAL ENTRY IS NOT LARGE ENOUGH. */ /* N - MUST BE SET TO THE ORDER OF THE MATRIX. IT IS NOT ALTERED. */ /* NZ - MUST BE SET TO THE NUMBER OF NON-ZEROS IN UPPER TRIANGLE OF */ /* PERMUTED MATRIX. NOT ALTERED BY MA27O/OD. */ /* A - MUST BE SET ON INPUT TO MATRIX HELD BY ROWS REORDERED BY */ /* PERMUTATION FROM MA27A/AD IN A(LA-NZ+I),I=1,NZ. ON */ /* EXIT FROM MA27O/OD, THE FACTORS OF U AND D ARE HELD IN */ /* POSITIONS 1 TO POSFAC-1. */ /* LA - LENGTH OF ARRAY A. A VALUE FOR LA */ /* SUFFICIENT FOR DEFINITE SYSTEMS */ /* WILL HAVE BEEN PROVIDED BY MA27A/AD. NOT ALTERED BY MA27O/OD. */ /* IW - MUST BE SET SO THAT,ON INPUT, IW(LIW-NZ+I),I=1,NZ */ /* HOLDS THE COLUMN INDEX OF THE ENTRY IN A(LA-NZ+I). ON EXIT, */ /* IW HOLDS INTEGER INDEXING INFORMATION ON THE FACTORS. */ /* THE ABSOLUTE VALUE OF THE FIRST ENTRY IN IW WILL BE SET TO */ /* THE NUMBER OF BLOCK PIVOTS ACTUALLY USED. THIS MAY BE */ /* DIFFERENT FROM NSTEPS SINCE NUMERICAL CONSIDERATIONS */ /* MAY PREVENT US CHOOSING A PIVOT AT EACH STAGE. IF THIS ENTRY */ /* IN IW IS NEGATIVE, THEN AT LEAST ONE TWO BY TWO */ /* PIVOT HAS BEEN USED DURING THE DECOMPOSITION. */ /* INTEGER INFORMATION ON EACH BLOCK PIVOT ROW FOLLOWS. FOR */ /* EACH BLOCK PIVOT ROW THE COLUMN INDICES ARE PRECEDED BY A */ /* COUNT OF THE NUMBER OF ROWS AND COLUMNS IN THE BLOCK PIVOT */ /* WHERE, IF ONLY ONE ROW IS PRESENT, ONLY THE NUMBER OF */ /* COLUMNS TOGETHER WITH A NEGATIVE FLAG IS HELD. THE FIRST */ /* COLUMN INDEX FOR A TWO BY TWO PIVOT IS FLAGGED NEGATIVE. */ /* LIW - LENGTH OF ARRAY IW. A VALUE FOR LIW SUFFICIENT FOR */ /* DEFINITE SYSTEMS */ /* WILL HAVE BEEN PROVIDED BY MA27A/AD. NOT ALTERED BY MA27O/OD */ /* PERM - PERM(I) MUST BE SET TO POSITION OF ROW/COLUMN I IN THE */ /* TENTATIVE PIVOT ORDER GENERATED BY MA27A/AD. */ /* IT IS NOT ALTERED BY MA27O/OD. */ /* NSTK - MUST BE LEFT UNCHANGED SINCE OUTPUT FROM MA27A/AD. NSTK(I) */ /* GIVES THE NUMBER OF GENERATED STACK ELEMENTS ASSEMBLED AT */ /* STAGE I. IT IS NOT ALTERED BY MA27O/OD. */ /* NSTEPS - LENGTH OF ARRAYS NSTK AND NELIM. VALUE IS GIVEN ON OUTPUT */ /* FROM MA27A/AD (WILL NEVER EXCEED N). IT IS NOT ALTERED BY */ /* MA27O/OD. */ /* MAXFRT - NEED NOT BE SET ON INPUT. ON OUTPUT */ /* MAXFRT WILL BE SET TO THE MAXIMUM FRONT SIZE ENCOUNTERED */ /* DURING THE DECOMPOSITION. */ /* NELIM - MUST BE UNCHANGED SINCE OUTPUT FROM MA27A/AD. NELIM(I) */ /* GIVES THE NUMBER OF ORIGINAL ROWS ASSEMBLED AT STAGE I. */ /* IT IS NOT ALTERED BY MA27O/OD. */ /* IW2 - INTEGER ARRAY OF LENGTH N. USED AS WORKSPACE BY MA27O/OD. */ /* ALTHOUGH WE COULD HAVE USED A SHORT WORD INTEGER IN THE IBM */ /* VERSION, WE HAVE NOT DONE SO BECAUSE THERE IS A SPARE */ /* FULL INTEGER ARRAY (USED IN THE SORT BEFORE MA27O/OD) */ /* AVAILABLE WHEN MA27O/OD IS CALLED FROM MA27B/BD. */ /* ICNTL is an INTEGER array of length 30, see MA27A/AD. */ /* CNTL is a REAL array of length 5, see MA27A/AD. */ /* INFO is an INTEGER array of length 20, see MA27A/AD. */ /* INFO(1) - INTEGER VARIABLE. DIAGNOSTIC FLAG. A ZERO VALUE ON EXIT */ /* INDICATES SUCCESS. POSSIBLE NEGATIVE VALUES ARE ... */ /* -3 INSUFFICIENT STORAGE FOR IW. */ /* -4 INSUFFICIENT STORAGE FOR A. */ /* -5 ZERO PIVOT FOUND IN FACTORIZATION OF DEFINITE MATRIX. */ /* .. Parameters .. */ /* .. */ /* .. Scalar Arguments .. */ /* .. */ /* .. Array Arguments .. */ /* .. */ /* .. Local Scalars .. */ /* .. */ /* .. External Subroutines .. */ /* .. */ /* .. Intrinsic Functions .. */ /* .. */ /* .. Statement Functions .. */ /* .. */ /* .. Statement Function definitions .. */ /* THE FOLLOWING ARITHMETIC FUNCTION GIVES THE DISPLACEMENT FROM */ /* THE START OF THE ASSEMBLED MATRIX(OF ORDER IX) OF THE DIAGONAL */ /* ENTRY IN ITS ROW IY. */ /* .. */ /* .. Executable Statements .. */ /* INITIALIZATION. */ /* NBLK IS THE NUMBER OF BLOCK PIVOTS USED. */ /* Parameter adjustments */ --iw2; --perm; --a; --iw; --nelim; --nstk; --icntl; --cntl; --info; /* Function Body */ nblk = 0; ntwo = 0; neig = 0; ncmpbi = 0; ncmpbr = 0; *maxfrt = 0; nrlbdu = 0; nirbdu = 0; /* A PRIVATE VARIABLE UU IS SET TO CNTL(1), SO THAT CNTL(1) WILL REMAIN */ /* UNALTERED. */ uu = dmin(cntl[1],.5f); uu = dmax(uu,-.5f); i__1 = *n; for (i__ = 1; i__ <= i__1; ++i__) { iw2[i__] = 0; /* L10: */ } /* IWPOS IS POINTER TO FIRST FREE POSITION FOR FACTORS IN IW. */ /* POSFAC IS POINTER FOR FACTORS IN A. AT EACH PASS THROUGH THE */ /* MAJOR LOOP POSFAC INITIALLY POINTS TO THE FIRST FREE LOCATION */ /* IN A AND THEN IS SET TO THE POSITION OF THE CURRENT PIVOT IN A. */ /* ISTK IS POINTER TO TOP OF STACK IN IW. */ /* ISTK2 IS POINTER TO BOTTOM OF STACK IN IW (NEEDED BY COMPRESS). */ /* ASTK IS POINTER TO TOP OF STACK IN A. */ /* ASTK2 IS POINTER TO BOTTOM OF STACK IN A (NEEDED BY COMPRESS). */ /* IINPUT IS POINTER TO CURRENT POSITION IN ORIGINAL ROWS IN IW. */ /* AINPUT IS POINTER TO CURRENT POSITION IN ORIGINAL ROWS IN A. */ /* AZERO IS POINTER TO LAST POSITION ZEROED IN A. */ /* NTOTPV IS THE TOTAL NUMBER OF PIVOTS SELECTED. THIS IS USED */ /* TO DETERMINE WHETHER THE MATRIX IS SINGULAR. */ iwpos = 2; posfac = 1; istk = *liw - *nz + 1; istk2 = istk - 1; astk = *la - *nz + 1; astk2 = astk - 1; iinput = istk; ainput = astk; azero = 0; ntotpv = 0; /* NUMASS IS THE ACCUMULATED NUMBER OF ROWS ASSEMBLED SO FAR. */ numass = 0; /* EACH PASS THROUGH THIS MAIN LOOP PERFORMS ALL THE OPERATIONS */ /* ASSOCIATED WITH ONE SET OF ASSEMBLY/ELIMINATIONS. */ i__1 = *nsteps; for (iass = 1; iass <= i__1; ++iass) { /* NASS WILL BE SET TO THE NUMBER OF FULLY ASSEMBLED VARIABLES IN */ /* CURRENT NEWLY CREATED ELEMENT. */ nass = nelim[iass]; /* NEWEL IS A POINTER INTO IW TO CONTROL OUTPUT OF INTEGER INFORMATION */ /* FOR NEWLY CREATED ELEMENT. */ newel = iwpos + 1; /* SYMBOLICALLY ASSEMBLE INCOMING ROWS AND GENERATED STACK ELEMENTS */ /* ORDERING THE RESULTANT ELEMENT ACCORDING TO PERMUTATION PERM. WE */ /* ASSEMBLE THE STACK ELEMENTS FIRST BECAUSE THESE WILL ALREADY BE */ /* ORDERED. */ /* SET HEADER POINTER FOR MERGE OF INDEX LISTS. */ jfirst = *n + 1; /* INITIALIZE NUMBER OF VARIABLES IN CURRENT FRONT. */ nfront = 0; numstk = nstk[iass]; ltopst = 1; lnass = 0; /* JUMP IF NO STACK ELEMENTS ARE BEING ASSEMBLED AT THIS STAGE. */ if (numstk == 0) { goto L80; } j2 = istk - 1; lnass = nass; ltopst = (iw[istk] + 1) * iw[istk] / 2; i__2 = numstk; for (iell = 1; iell <= i__2; ++iell) { /* ASSEMBLE ELEMENT IELL PLACING */ /* THE INDICES INTO A LINKED LIST IN IW2 ORDERED */ /* ACCORDING TO PERM. */ jnext = jfirst; jlast = *n + 1; j1 = j2 + 2; j2 = j1 - 1 + iw[j1 - 1]; /* RUN THROUGH INDEX LIST OF STACK ELEMENT IELL. */ i__3 = j2; for (jj = j1; jj <= i__3; ++jj) { j = iw[jj]; /* JUMP IF ALREADY ASSEMBLED */ if (iw2[j] > 0) { goto L60; } jnew = perm[j]; /* IF VARIABLE WAS PREVIOUSLY FULLY SUMMED BUT WAS NOT PIVOTED ON */ /* EARLIER BECAUSE OF NUMERICAL TEST, INCREMENT NUMBER OF FULLY */ /* SUMMED ROWS/COLUMNS IN FRONT. */ if (jnew <= numass) { ++nass; } /* FIND POSITION IN LINKED LIST FOR NEW VARIABLE. NOTE THAT WE START */ /* FROM WHERE WE LEFT OFF AFTER ASSEMBLY OF PREVIOUS VARIABLE. */ i__4 = *n; for (idummy = 1; idummy <= i__4; ++idummy) { if (jnext == *n + 1) { goto L30; } if (perm[jnext] > jnew) { goto L30; } jlast = jnext; jnext = iw2[jlast]; /* L20: */ } L30: if (jlast != *n + 1) { goto L40; } jfirst = j; goto L50; L40: iw2[jlast] = j; L50: iw2[j] = jnext; jlast = j; /* INCREMENT NUMBER OF VARIABLES IN THE FRONT. */ ++nfront; L60: ; } /* L70: */ } lnass = nass - lnass; /* NOW INCORPORATE ORIGINAL ROWS. NOTE THAT THE COLUMNS IN THESE ROWS */ /* NEED NOT BE IN ORDER. WE ALSO PERFORM */ /* A SWOP SO THAT THE DIAGONAL ENTRY IS THE FIRST IN ITS */ /* ROW. THIS ALLOWS US TO AVOID STORING THE INVERSE OF ARRAY PERM. */ L80: numorg = nelim[iass]; j1 = iinput; i__2 = numorg; for (iorg = 1; iorg <= i__2; ++iorg) { j = -iw[j1]; i__3 = *liw; for (idummy = 1; idummy <= i__3; ++idummy) { jnew = perm[j]; /* JUMP IF VARIABLE ALREADY INCLUDED. */ if (iw2[j] > 0) { goto L130; } /* HERE WE MUST ALWAYS START OUR SEARCH AT THE BEGINNING. */ jlast = *n + 1; jnext = jfirst; i__4 = *n; for (jdummy = 1; jdummy <= i__4; ++jdummy) { if (jnext == *n + 1) { goto L100; } if (perm[jnext] > jnew) { goto L100; } jlast = jnext; jnext = iw2[jlast]; /* L90: */ } L100: if (jlast != *n + 1) { goto L110; } jfirst = j; goto L120; L110: iw2[jlast] = j; L120: iw2[j] = jnext; /* INCREMENT NUMBER OF VARIABLES IN FRONT. */ ++nfront; L130: ++j1; if (j1 > *liw) { goto L150; } j = iw[j1]; if (j < 0) { goto L150; } /* L140: */ } L150: ; } /* NOW RUN THROUGH LINKED LIST IW2 PUTTING INDICES OF VARIABLES IN NEW */ /* ELEMENT INTO IW AND SETTING IW2 ENTRY TO POINT TO THE RELATIVE */ /* POSITION OF THE VARIABLE IN THE NEW ELEMENT. */ if (newel + nfront < istk) { goto L160; } /* COMPRESS IW. */ ma27p_(&a[1], &iw[1], &istk, &istk2, &iinput, &c__2, &ncmpbr, &ncmpbi) ; if (newel + nfront < istk) { goto L160; } info[2] = *liw + 1 + newel + nfront - istk; goto L770; L160: j = jfirst; i__2 = nfront; for (ifr = 1; ifr <= i__2; ++ifr) { ++newel; iw[newel] = j; jnext = iw2[j]; iw2[j] = newel - (iwpos + 1); j = jnext; /* L170: */ } /* ASSEMBLE REALS INTO FRONTAL MATRIX. */ *maxfrt = max(*maxfrt,nfront); iw[iwpos] = nfront; /* FIRST ZERO OUT FRONTAL MATRIX AS APPROPRIATE FIRST CHECKING TO SEE */ /* IF THERE IS SUFFICIENT SPACE. */ laell = (nfront + 1) * nfront / 2; apos2 = posfac + laell - 1; if (numstk != 0) { lnass = lnass * ((nfront << 1) - lnass + 1) / 2; } if (posfac + lnass - 1 >= astk) { goto L180; } if (apos2 < astk + ltopst - 1) { goto L190; } /* COMPRESS A. */ L180: ma27p_(&a[1], &iw[1], &astk, &astk2, &ainput, &c__1, &ncmpbr, &ncmpbi) ; if (posfac + lnass - 1 >= astk) { goto L780; } if (apos2 >= astk + ltopst - 1) { goto L780; } L190: if (apos2 <= azero) { goto L220; } apos = azero + 1; /* Computing MIN */ i__2 = apos2, i__3 = astk - 1; lapos2 = min(i__2,i__3); if (lapos2 < apos) { goto L210; } i__2 = lapos2; for (k = apos; k <= i__2; ++k) { a[k] = 0.f; /* L200: */ } L210: azero = apos2; /* JUMP IF THERE ARE NO STACK ELEMENTS TO ASSEMBLE. */ L220: if (numstk == 0) { goto L260; } /* PLACE REALS CORRESPONDING TO STACK ELEMENTS IN CORRECT POSITIONS IN A. */ i__2 = numstk; for (iell = 1; iell <= i__2; ++iell) { j1 = istk + 1; j2 = istk + iw[istk]; i__3 = j2; for (jj = j1; jj <= i__3; ++jj) { irow = iw[jj]; irow = iw2[irow]; apos = posfac + (irow - 1) * (2 * nfront - irow + 2) / 2; i__4 = j2; for (jjj = jj; jjj <= i__4; ++jjj) { j = iw[jjj]; apos2 = apos + iw2[j] - irow; a[apos2] += a[astk]; a[astk] = 0.f; ++astk; /* L230: */ } /* L240: */ } /* INCREMENT STACK POINTER. */ istk = j2 + 1; /* L250: */ } /* INCORPORATE REALS FROM ORIGINAL ROWS. */ L260: i__2 = numorg; for (iorg = 1; iorg <= i__2; ++iorg) { j = -iw[iinput]; /* WE CAN DO THIS BECAUSE THE DIAGONAL IS NOW THE FIRST ENTRY. */ irow = iw2[j]; apos = posfac + (irow - 1) * (2 * nfront - irow + 2) / 2; /* THE FOLLOWING LOOP GOES FROM 1 TO NZ BECAUSE THERE MAY BE DUPLICATES. */ i__3 = *nz; for (idummy = 1; idummy <= i__3; ++idummy) { apos2 = apos + iw2[j] - irow; a[apos2] += a[ainput]; ++ainput; ++iinput; if (iinput > *liw) { goto L280; } j = iw[iinput]; if (j < 0) { goto L280; } /* L270: */ } L280: ; } /* RESET IW2 AND NUMASS. */ numass += numorg; j1 = iwpos + 2; j2 = iwpos + nfront + 1; i__2 = j2; for (k = j1; k <= i__2; ++k) { j = iw[k]; iw2[j] = 0; /* L290: */ } /* PERFORM PIVOTING ON ASSEMBLED ELEMENT. */ /* NPIV IS THE NUMBER OF PIVOTS SO FAR SELECTED. */ /* LNPIV IS THE NUMBER OF PIVOTS SELECTED AFTER THE LAST PASS THROUGH */ /* THE THE FOLLOWING LOOP. */ lnpiv = -1; npiv = 0; i__2 = nass; for (kdummy = 1; kdummy <= i__2; ++kdummy) { if (npiv == nass) { goto L660; } if (npiv == lnpiv) { goto L660; } lnpiv = npiv; npivp1 = npiv + 1; /* JPIV IS USED AS A FLAG TO INDICATE WHEN 2 BY 2 PIVOTING HAS OCCURRED */ /* SO THAT IPIV IS INCREMENTED CORRECTLY. */ jpiv = 1; /* NASS IS MAXIMUM POSSIBLE NUMBER OF PIVOTS. */ /* WE EITHER TAKE THE DIAGONAL ENTRY OR THE 2 BY 2 PIVOT WITH THE */ /* LARGEST OFF-DIAGONAL AT EACH STAGE. */ /* EACH PASS THROUGH THIS LOOP TRIES TO CHOOSE ONE PIVOT. */ i__3 = nass; for (ipiv = npivp1; ipiv <= i__3; ++ipiv) { --jpiv; /* JUMP IF WE HAVE JUST PROCESSED A 2 BY 2 PIVOT. */ if (jpiv == 1) { goto L640; } i__4 = nfront - npiv; i__5 = ipiv - npiv; apos = posfac + (i__5 - 1) * (2 * i__4 - i__5 + 2) / 2; /* IF THE USER HAS INDICATED THAT THE MATRIX IS DEFINITE, WE */ /* DO NOT NEED TO TEST FOR STABILITY BUT WE DO CHECK TO SEE IF THE */ /* PIVOT IS NON-ZERO OR HAS CHANGED SIGN. */ /* IF IT IS ZERO, WE EXIT WITH AN ERROR. IF IT HAS CHANGED SIGN */ /* AND U WAS SET NEGATIVE, THEN WE AGAIN EXIT IMMEDIATELY. IF THE */ /* PIVOT CHANGES SIGN AND U WAS ZERO, WE CONTINUE WITH THE */ /* FACTORIZATION BUT PRINT A WARNING MESSAGE ON UNIT ICNTL(2). */ /* ISNPIV HOLDS A FLAG FOR THE SIGN OF THE PIVOTS TO DATE SO THAT */ /* A SIGN CHANGE WHEN DECOMPOSING AN ALLEGEDLY DEFINITE MATRIX CAN */ /* BE DETECTED. */ if (uu > 0.f) { goto L320; } if ((r__1 = a[apos], dabs(r__1)) <= cntl[3]) { goto L790; } /* JUMP IF THIS IS NOT THE FIRST PIVOT TO BE SELECTED. */ if (ntotpv > 0) { goto L300; } /* SET ISNPIV. */ if (a[apos] > 0.f) { isnpiv = 1; } if (a[apos] < 0.f) { isnpiv = -1; } L300: if (a[apos] > 0.f && isnpiv == 1) { goto L560; } if (a[apos] < 0.f && isnpiv == -1) { goto L560; } if (info[1] != 2) { info[2] = 0; } ++info[2]; info[1] = 2; i__ = ntotpv + 1; if (icntl[2] > 0 && info[2] <= 10) { io___280.ciunit = icntl[2]; s_wsfe(&io___280); do_fio(&c__1, (char *)&info[1], (ftnlen)sizeof(integer)); do_fio(&c__1, (char *)&i__, (ftnlen)sizeof(integer)); e_wsfe(); } isnpiv = -isnpiv; if (uu == 0.f) { goto L560; } goto L800; L320: amax = 0.f; tmax = amax; /* FIND LARGEST ENTRY TO RIGHT OF DIAGONAL IN ROW OF PROSPECTIVE PIVOT */ /* IN THE FULLY-SUMMED PART. ALSO RECORD COLUMN OF THIS LARGEST */ /* ENTRY. */ j1 = apos + 1; j2 = apos + nass - ipiv; if (j2 < j1) { goto L340; } i__4 = j2; for (jj = j1; jj <= i__4; ++jj) { if ((r__1 = a[jj], dabs(r__1)) <= amax) { goto L330; } jmax = ipiv + jj - j1 + 1; amax = (r__1 = a[jj], dabs(r__1)); L330: ; } /* DO SAME AS ABOVE FOR NON-FULLY-SUMMED PART ONLY HERE WE DO NOT NEED */ /* TO RECORD COLUMN SO LOOP IS SIMPLER. */ L340: j1 = j2 + 1; j2 = apos + nfront - ipiv; if (j2 < j1) { goto L360; } i__4 = j2; for (jj = j1; jj <= i__4; ++jj) { /* Computing MAX */ r__2 = (r__1 = a[jj], dabs(r__1)); tmax = dmax(r__2,tmax); /* L350: */ } /* NOW CALCULATE LARGEST ENTRY IN OTHER PART OF ROW. */ L360: rmax = dmax(tmax,amax); apos1 = apos; kk = nfront - ipiv; lt = ipiv - (npiv + 1); if (lt == 0) { goto L380; } i__4 = lt; for (k = 1; k <= i__4; ++k) { ++kk; apos1 -= kk; /* Computing MAX */ r__2 = rmax, r__3 = (r__1 = a[apos1], dabs(r__1)); rmax = dmax(r__2,r__3); /* L370: */ } /* JUMP IF STABILITY TEST SATISFIED. */ L380: /* Computing MAX */ r__2 = cntl[3], r__3 = uu * rmax; if ((r__1 = a[apos], dabs(r__1)) > dmax(r__2,r__3)) { goto L450; } /* CHECK BLOCK PIVOT OF ORDER 2 FOR STABILITY. */ if (dabs(amax) <= cntl[3]) { goto L640; } i__4 = nfront - npiv; i__5 = jmax - npiv; apos2 = posfac + (i__5 - 1) * (2 * i__4 - i__5 + 2) / 2; detpiv = a[apos] * a[apos2] - amax * amax; thresh = dabs(detpiv); /* SET THRESH TO U TIMES THE RECIPROCAL OF THE MAX-NORM OF THE INVERSE */ /* OF THE PROSPECTIVE BLOCK. */ /* Computing MAX */ r__3 = (r__1 = a[apos], dabs(r__1)) + amax, r__4 = (r__2 = a[ apos2], dabs(r__2)) + amax; thresh /= uu * dmax(r__3,r__4); /* CHECK 2 BY 2 PIVOT FOR STABILITY. */ /* FIRST CHECK AGAINST ROW IPIV. */ if (thresh <= rmax) { goto L640; } /* FIND LARGEST ENTRY IN ROW JMAX. */ /* FIND MAXIMUM TO THE RIGHT OF THE DIAGONAL. */ rmax = 0.f; j1 = apos2 + 1; j2 = apos2 + nfront - jmax; if (j2 < j1) { goto L400; } i__4 = j2; for (jj = j1; jj <= i__4; ++jj) { /* Computing MAX */ r__2 = rmax, r__3 = (r__1 = a[jj], dabs(r__1)); rmax = dmax(r__2,r__3); /* L390: */ } /* NOW CHECK TO THE LEFT OF THE DIAGONAL. */ /* WE USE TWO LOOPS TO AVOID TESTING FOR ROW IPIV INSIDE THE LOOP. */ L400: kk = nfront - jmax + 1; apos3 = apos2; jmxmip = jmax - ipiv - 1; if (jmxmip == 0) { goto L420; } i__4 = jmxmip; for (k = 1; k <= i__4; ++k) { apos2 -= kk; ++kk; /* Computing MAX */ r__2 = rmax, r__3 = (r__1 = a[apos2], dabs(r__1)); rmax = dmax(r__2,r__3); /* L410: */ } L420: ipmnp = ipiv - npiv - 1; if (ipmnp == 0) { goto L440; } apos2 -= kk; ++kk; i__4 = ipmnp; for (k = 1; k <= i__4; ++k) { apos2 -= kk; ++kk; /* Computing MAX */ r__2 = rmax, r__3 = (r__1 = a[apos2], dabs(r__1)); rmax = dmax(r__2,r__3); /* L430: */ } L440: if (thresh <= rmax) { goto L640; } pivsiz = 2; goto L460; L450: pivsiz = 1; L460: irow = ipiv - npiv; /* PIVOT HAS BEEN CHOSEN. IF BLOCK PIVOT OF ORDER 2, PIVSIZ IS EQUAL TO */ /* TWO OTHERWISE PIVSIZ EQUALS ONE.. */ /* THE FOLLOWING LOOP MOVES THE PIVOT BLOCK TO THE TOP LEFT HAND CORNER */ /* OF THE FRONTAL MATRIX. */ i__4 = pivsiz; for (krow = 1; krow <= i__4; ++krow) { /* WE JUMP IF SWOP IS NOT NECESSARY. */ if (irow == 1) { goto L530; } j1 = posfac + irow; j2 = posfac + nfront - (npiv + 1); if (j2 < j1) { goto L480; } apos2 = apos + 1; /* SWOP PORTION OF ROWS WHOSE COLUMN INDICES ARE GREATER THAN LATER ROW. */ i__5 = j2; for (jj = j1; jj <= i__5; ++jj) { swop = a[apos2]; a[apos2] = a[jj]; a[jj] = swop; ++apos2; /* L470: */ } L480: j1 = posfac + 1; j2 = posfac + irow - 2; apos2 = apos; kk = nfront - (irow + npiv); if (j2 < j1) { goto L500; } /* SWOP PORTION OF ROWS/COLUMNS WHOSE INDICES LIE BETWEEN THE TWO ROWS. */ i__5 = j2; for (jjj = j1; jjj <= i__5; ++jjj) { jj = j2 - jjj + j1; ++kk; apos2 -= kk; swop = a[apos2]; a[apos2] = a[jj]; a[jj] = swop; /* L490: */ } L500: if (npiv == 0) { goto L520; } apos1 = posfac; ++kk; apos2 -= kk; /* SWOP PORTION OF COLUMNS WHOSE INDICES ARE LESS THAN EARLIER ROW. */ i__5 = npiv; for (jj = 1; jj <= i__5; ++jj) { ++kk; apos1 -= kk; apos2 -= kk; swop = a[apos2]; a[apos2] = a[apos1]; a[apos1] = swop; /* L510: */ } /* SWOP DIAGONALS AND INTEGER INDEXING INFORMATION */ L520: swop = a[apos]; a[apos] = a[posfac]; a[posfac] = swop; ipos = iwpos + npiv + 2; iexch = iwpos + irow + npiv + 1; iswop = iw[ipos]; iw[ipos] = iw[iexch]; iw[iexch] = iswop; L530: if (pivsiz == 1) { goto L550; } /* SET VARIABLES FOR THE SWOP OF SECOND ROW OF BLOCK PIVOT. */ if (krow == 2) { goto L540; } irow = jmax - (npiv + 1); jpos = posfac; posfac = posfac + nfront - npiv; ++npiv; apos = apos3; goto L550; /* RESET VARIABLES PREVIOUSLY SET FOR SECOND PASS. */ L540: --npiv; posfac = jpos; L550: ; } if (pivsiz == 2) { goto L600; } /* PERFORM THE ELIMINATION USING ENTRY (IPIV,IPIV) AS PIVOT. */ /* WE STORE U AND DINVERSE. */ L560: a[posfac] = 1.f / a[posfac]; if (a[posfac] < 0.f) { ++neig; } j1 = posfac + 1; j2 = posfac + nfront - (npiv + 1); if (j2 < j1) { goto L590; } ibeg = j2 + 1; i__4 = j2; for (jj = j1; jj <= i__4; ++jj) { amult = -a[jj] * a[posfac]; iend = ibeg + nfront - (npiv + jj - j1 + 2); /* THE FOLLOWING SPECIAL COMMENT FORCES VECTORIZATION ON THE CRAY-1. */ /* DIR$ IVDEP */ i__5 = iend; for (irow = ibeg; irow <= i__5; ++irow) { jcol = jj + irow - ibeg; a[irow] += amult * a[jcol]; /* L570: */ } ibeg = iend + 1; a[jj] = amult; /* L580: */ } L590: ++npiv; ++ntotpv; jpiv = 1; posfac = posfac + nfront - npiv + 1; goto L640; /* PERFORM ELIMINATION USING BLOCK PIVOT OF ORDER TWO. */ /* REPLACE BLOCK PIVOT BY ITS INVERSE. */ /* SET FLAG TO INDICATE USE OF 2 BY 2 PIVOT IN IW. */ L600: ipos = iwpos + npiv + 2; ++ntwo; iw[ipos] = -iw[ipos]; pospv1 = posfac; pospv2 = posfac + nfront - npiv; swop = a[pospv2]; if (detpiv < 0.f) { ++neig; } if (detpiv > 0.f && swop < 0.f) { neig += 2; } a[pospv2] = a[pospv1] / detpiv; a[pospv1] = swop / detpiv; a[pospv1 + 1] = -a[pospv1 + 1] / detpiv; j1 = pospv1 + 2; j2 = pospv1 + nfront - (npiv + 1); if (j2 < j1) { goto L630; } jj1 = pospv2; ibeg = pospv2 + nfront - (npiv + 1); i__4 = j2; for (jj = j1; jj <= i__4; ++jj) { ++jj1; amult1 = -(a[pospv1] * a[jj] + a[pospv1 + 1] * a[jj1]); amult2 = -(a[pospv1 + 1] * a[jj] + a[pospv2] * a[jj1]); iend = ibeg + nfront - (npiv + jj - j1 + 3); /* THE FOLLOWING SPECIAL COMMENT FORCES VECTORIZATION ON THE CRAY-1. */ /* DIR$ IVDEP */ i__5 = iend; for (irow = ibeg; irow <= i__5; ++irow) { k1 = jj + irow - ibeg; k2 = jj1 + irow - ibeg; a[irow] = a[irow] + amult1 * a[k1] + amult2 * a[k2]; /* L610: */ } ibeg = iend + 1; a[jj] = amult1; a[jj1] = amult2; /* L620: */ } L630: npiv += 2; ntotpv += 2; jpiv = 2; posfac = pospv2 + nfront - npiv + 1; L640: ; } /* L650: */ } /* END OF MAIN ELIMINATION LOOP. */ L660: if (npiv != 0) { ++nblk; } ioldps = iwpos; iwpos = iwpos + nfront + 2; if (npiv == 0) { goto L690; } if (npiv > 1) { goto L680; } iw[ioldps] = -iw[ioldps]; i__2 = nfront; for (k = 1; k <= i__2; ++k) { j1 = ioldps + k; iw[j1] = iw[j1 + 1]; /* L670: */ } --iwpos; goto L690; L680: iw[ioldps + 1] = npiv; /* COPY REMAINDER OF ELEMENT TO TOP OF STACK */ L690: liell = nfront - npiv; if (liell == 0 || iass == *nsteps) { goto L750; } if (iwpos + liell < istk) { goto L700; } ma27p_(&a[1], &iw[1], &istk, &istk2, &iinput, &c__2, &ncmpbr, &ncmpbi) ; L700: istk = istk - liell - 1; iw[istk] = liell; j1 = istk; kk = iwpos - liell - 1; /* THE FOLLOWING SPECIAL COMMENT FORCES VECTORIZATION ON THE CRAY-1. */ /* DIR$ IVDEP */ i__2 = liell; for (k = 1; k <= i__2; ++k) { ++j1; ++kk; iw[j1] = iw[kk]; /* L710: */ } /* WE COPY IN REVERSE DIRECTION TO AVOID OVERWRITE PROBLEMS. */ laell = (liell + 1) * liell / 2; kk = posfac + laell; if (kk != astk) { goto L720; } astk -= laell; goto L740; /* THE MOVE AND ZEROING OF ARRAY A IS PERFORMED WITH TWO LOOPS SO */ /* THAT THE CRAY-1 WILL VECTORIZE THEM SAFELY. */ L720: kmax = kk - 1; /* THE FOLLOWING SPECIAL COMMENT FORCES VECTORIZATION ON THE CRAY-1. */ /* DIR$ IVDEP */ i__2 = laell; for (k = 1; k <= i__2; ++k) { --kk; --astk; a[astk] = a[kk]; /* L730: */ } /* Computing MIN */ i__2 = kmax, i__3 = astk - 1; kmax = min(i__2,i__3); i__2 = kmax; for (k = kk; k <= i__2; ++k) { a[k] = 0.f; /* L735: */ } L740: /* Computing MIN */ i__2 = azero, i__3 = astk - 1; azero = min(i__2,i__3); L750: if (npiv == 0) { iwpos = ioldps; } /* L760: */ } /* END OF LOOP ON TREE NODES. */ iw[1] = nblk; if (ntwo > 0) { iw[1] = -nblk; } nrlbdu = posfac - 1; nirbdu = iwpos - 1; if (ntotpv == *n) { goto L810; } info[1] = 3; info[2] = ntotpv; goto L810; /* **** ERROR RETURNS **** */ L770: info[1] = -3; goto L810; L780: info[1] = -4; /* Computing MAX */ i__1 = posfac + lnass, i__2 = apos2 - ltopst + 2; info[2] = *la + max(i__1,i__2) - astk; goto L810; L790: info[1] = -5; info[2] = ntotpv + 1; goto L810; L800: info[1] = -6; info[2] = ntotpv + 1; L810: info[9] = nrlbdu; info[10] = nirbdu; info[12] = ncmpbr; info[13] = ncmpbi; info[14] = ntwo; info[15] = neig; return 0; } /* ma27o_ */ /* Subroutine */ int ma27p_(real *a, integer *iw, integer *j1, integer *j2, integer *itop, integer *ireal, integer *ncmpbr, integer *ncmpbi) { /* System generated locals */ integer i__1; /* Local variables */ static integer jj, jjj, ipos; /* THIS SUBROUTINE PERFORMS A VERY SIMPLE COMPRESS (BLOCK MOVE). */ /* ENTRIES J1 TO J2 (INCL.) IN A OR IW AS APPROPRIATE ARE MOVED TO */ /* OCCUPY THE POSITIONS IMMEDIATELY PRIOR TO POSITION ITOP. */ /* A/IW HOLD THE ARRAY BEING COMPRESSED. */ /* J1/J2 DEFINE THE ENTRIES BEING MOVED. */ /* ITOP DEFINES THE POSITION IMMEDIATELY AFTER THE POSITIONS TO WHICH */ /* J1 TO J2 ARE MOVED. */ /* IREAL MUST BE SET BY THE USER TO 2 IF THE MOVE IS ON ARRAY IW, */ /* ANY OTHER VALUE WILL PERFORM THE MOVE ON A. */ /* NCMPBR and NCMPBI, see INFO(12) and INFO(13) in MA27A/AD (ACCUMULATE */ /* THE NUMBER OF COMPRESSES OF THE REALS AND INTEGERS PERFORMED BY */ /* MA27B/BD. */ /* .. Scalar Arguments .. */ /* .. */ /* .. Array Arguments .. */ /* .. */ /* .. Local Scalars .. */ /* .. */ /* .. Executable Statements .. */ /* Parameter adjustments */ --iw; --a; /* Function Body */ ipos = *itop - 1; if (*j2 == ipos) { goto L50; } if (*ireal == 2) { goto L20; } ++(*ncmpbr); if (*j1 > *j2) { goto L40; } i__1 = *j2; for (jjj = *j1; jjj <= i__1; ++jjj) { jj = *j2 - jjj + *j1; a[ipos] = a[jj]; --ipos; /* L10: */ } goto L40; L20: ++(*ncmpbi); if (*j1 > *j2) { goto L40; } i__1 = *j2; for (jjj = *j1; jjj <= i__1; ++jjj) { jj = *j2 - jjj + *j1; iw[ipos] = iw[jj]; --ipos; /* L30: */ } L40: *j2 = *itop - 1; *j1 = ipos + 1; L50: return 0; } /* ma27p_ */ /* Subroutine */ int ma27q_(integer *n, real *a, integer *la, integer *iw, integer *liw, real *w, integer *maxfnt, real *rhs, integer *iw2, integer *nblk, integer *latop, integer *icntl) { /* System generated locals */ integer i__1, i__2, i__3; /* Local variables */ static integer j, k, j1, j2, j3, k1, k2, k3; static real w1, w2; static integer jj, ifr, ist, iblk, apos, irhs, ilvl, ipiv, jpiv, ipos, npiv, irow, liell; /* THIS SUBROUTINE PERFORMS FORWARD ELIMINATION */ /* USING THE FACTOR U TRANSPOSE STORED IN A/IA AFTER MA27B/BD. */ /* N - MUST BE SET TO THE ORDER OF THE MATRIX. NOT ALTERED */ /* BY MA27Q/QD. */ /* A - MUST BE SET TO HOLD THE REAL VALUES */ /* CORRESPONDING TO THE FACTORS OF DINVERSE AND U. THIS MUST BE */ /* UNCHANGED SINCE THE PRECEDING CALL TO MA27B/BD. NOT ALTERED */ /* BY MA27Q/QD. */ /* LA - LENGTH OF ARRAY A. NOT ALTERED BY MA27Q/QD. */ /* IW - HOLDS THE INTEGER INDEXING */ /* INFORMATION FOR THE MATRIX FACTORS IN A. THIS MUST BE */ /* UNCHANGED SINCE THE PRECEDING CALL TO MA27B/BD. NOT ALTERED */ /* BY MA27Q/QD. */ /* LIW - LENGTH OF ARRAY IW. NOT ALTERED BY MA27Q/QD. */ /* W - USED */ /* AS WORKSPACE BY MA27Q/QD TO HOLD THE COMPONENTS OF THE RIGHT */ /* HAND SIDES CORRESPONDING TO CURRENT BLOCK PIVOTAL ROWS. */ /* MAXFNT - MUST BE SET TO THE LARGEST NUMBER OF */ /* VARIABLES IN ANY BLOCK PIVOT ROW. THIS VALUE WILL HAVE */ /* BEEN OUTPUT BY MA27B/BD. NOT ALTERED BY MA27Q/QD. */ /* RHS - ON INPUT, */ /* MUST BE SET TO HOLD THE RIGHT HAND SIDES FOR THE EQUATIONS */ /* WHICH THE USER DESIRES TO SOLVE. ON OUTPUT, RHS WILL HOLD */ /* THE MODIFIED VECTORS CORRESPONDING TO PERFORMING FORWARD */ /* ELIMINATION ON THE RIGHT HAND SIDES. */ /* IW2 - NEED NOT BE SET ON ENTRY. ON EXIT IW2(I) (I = 1,NBLK) */ /* WILL HOLD POINTERS TO THE */ /* BEGINNING OF EACH BLOCK PIVOT IN ARRAY IW. */ /* NBLK - NUMBER OF BLOCK PIVOT ROWS. NOT ALTERED BY MA27Q/QD. */ /* LATOP - NEED NOT BE SET ON ENTRY. ON EXIT, IT IS THE POSITION IN */ /* A OF THE LAST ENTRY IN THE FACTORS. IT MUST BE PASSED */ /* UNCHANGED TO MA27R/RD. */ /* ICNTL is an INTEGER array of length 30, see MA27A/AD. */ /* ICNTL(IFRLVL+I) I=1,20 IS USED TO CONTROL WHETHER DIRECT OR INDIRECT */ /* ACCESS IS USED BY MA27C/CD. INDIRECT ACCESS IS EMPLOYED */ /* IN FORWARD AND BACK SUBSTITUTION RESPECTIVELY IF THE SIZE OF */ /* A BLOCK IS LESS THAN ICNTL(IFRLVL+MIN(10,NPIV)) AND */ /* ICNTL(IFRLVL+10+MIN(10,NPIV)) RESPECTIVELY, WHERE NPIV IS THE */ /* NUMBER OF PIVOTS IN THE BLOCK. */ /* .. Scalar Arguments .. */ /* .. */ /* .. Array Arguments .. */ /* .. */ /* .. Local Scalars .. */ /* .. */ /* .. Intrinsic Functions .. */ /* .. */ /* .. Executable Statements .. */ /* APOS. RUNNING POINTER TO CURRENT PIVOT POSITION IN ARRAY A. */ /* IPOS. RUNNING POINTER TO BEGINNING OF BLOCK PIVOT ROW IN IW. */ /* Parameter adjustments */ --rhs; --a; --iw; --w; --iw2; --icntl; /* Function Body */ apos = 1; ipos = 1; j2 = 0; iblk = 0; npiv = 0; i__1 = *n; for (irow = 1; irow <= i__1; ++irow) { if (npiv > 0) { goto L90; } ++iblk; if (iblk > *nblk) { goto L150; } ipos = j2 + 1; /* SET UP POINTER IN PREPARATION FOR BACK SUBSTITUTION. */ iw2[iblk] = ipos; /* ABS(LIELL) IS NUMBER OF VARIABLES (COLUMNS) IN BLOCK PIVOT ROW. */ liell = -iw[ipos]; /* NPIV IS NUMBER OF PIVOTS (ROWS) IN BLOCK PIVOT. */ npiv = 1; if (liell > 0) { goto L10; } liell = -liell; ++ipos; npiv = iw[ipos]; L10: j1 = ipos + 1; j2 = ipos + liell; ilvl = min(npiv,10); if (liell < icntl[ilvl + 5]) { goto L90; } /* PERFORM OPERATIONS USING DIRECT ADDRESSING. */ /* LOAD APPROPRIATE COMPONENTS OF RIGHT HAND SIDES INTO ARRAY W. */ ifr = 0; i__2 = j2; for (jj = j1; jj <= i__2; ++jj) { j = (i__3 = iw[jj], abs(i__3)); ++ifr; w[ifr] = rhs[j]; /* L20: */ } /* JPIV IS USED AS A FLAG SO THAT IPIV IS INCREMENTED CORRECTLY AFTER */ /* THE USE OF A 2 BY 2 PIVOT. */ jpiv = 1; j3 = j1; /* PERFORM OPERATIONS. */ i__2 = npiv; for (ipiv = 1; ipiv <= i__2; ++ipiv) { --jpiv; if (jpiv == 1) { goto L70; } /* JUMP IF WE HAVE A 2 BY 2 PIVOT. */ if (iw[j3] < 0) { goto L40; } /* PERFORM FORWARD SUBSTITUTION USING 1 BY 1 PIVOT. */ jpiv = 1; ++j3; ++apos; ist = ipiv + 1; if (liell < ist) { goto L70; } w1 = w[ipiv]; k = apos; i__3 = liell; for (j = ist; j <= i__3; ++j) { w[j] += a[k] * w1; ++k; /* L30: */ } apos = apos + liell - ist + 1; goto L70; /* PERFORM OPERATIONS WITH 2 BY 2 PIVOT. */ L40: jpiv = 2; j3 += 2; apos += 2; ist = ipiv + 2; if (liell < ist) { goto L60; } w1 = w[ipiv]; w2 = w[ipiv + 1]; k1 = apos; k2 = apos + liell - ipiv; i__3 = liell; for (j = ist; j <= i__3; ++j) { w[j] = w[j] + w1 * a[k1] + w2 * a[k2]; ++k1; ++k2; /* L50: */ } L60: apos = apos + (liell - ist + 1 << 1) + 1; L70: ; } /* RELOAD W BACK INTO RHS. */ ifr = 0; i__2 = j2; for (jj = j1; jj <= i__2; ++jj) { j = (i__3 = iw[jj], abs(i__3)); ++ifr; rhs[j] = w[ifr]; /* L80: */ } npiv = 0; goto L140; /* PERFORM OPERATIONS USING INDIRECT ADDRESSING. */ /* JUMP IF WE HAVE A 2 BY 2 PIVOT. */ L90: if (iw[j1] < 0) { goto L110; } /* PERFORM FORWARD SUBSTITUTION USING 1 BY 1 PIVOT. */ --npiv; ++apos; ++j1; if (j1 > j2) { goto L140; } irhs = iw[j1 - 1]; w1 = rhs[irhs]; k = apos; i__2 = j2; for (j = j1; j <= i__2; ++j) { irhs = (i__3 = iw[j], abs(i__3)); rhs[irhs] += a[k] * w1; ++k; /* L100: */ } apos = apos + j2 - j1 + 1; goto L140; /* PERFORM OPERATIONS WITH 2 BY 2 PIVOT */ L110: npiv += -2; j1 += 2; apos += 2; if (j1 > j2) { goto L130; } irhs = -iw[j1 - 2]; w1 = rhs[irhs]; irhs = iw[j1 - 1]; w2 = rhs[irhs]; k1 = apos; k3 = apos + j2 - j1 + 2; i__2 = j2; for (j = j1; j <= i__2; ++j) { irhs = (i__3 = iw[j], abs(i__3)); rhs[irhs] = rhs[irhs] + w1 * a[k1] + w2 * a[k3]; ++k1; ++k3; /* L120: */ } L130: apos = apos + (j2 - j1 + 1 << 1) + 1; L140: ; } L150: *latop = apos - 1; return 0; } /* ma27q_ */ /* Subroutine */ int ma27r_(integer *n, real *a, integer *la, integer *iw, integer *liw, real *w, integer *maxfnt, real *rhs, integer *iw2, integer *nblk, integer *latop, integer *icntl) { /* System generated locals */ integer i__1, i__2, i__3; /* Local variables */ static integer j, k, j1, j2; static real w1, w2; static integer jj, jj1, jj2, ifr, ist, iblk, apos, irhs, ilvl, ipiv, jpiv, loop, ipos, jpos, npiv, apos2, i1rhs, i2rhs, liell, iirhs, iipiv; /* THIS SUBROUTINE PERFORMS BACKWARD ELIMINATION OPERATIONS */ /* USING THE FACTORS DINVERSE AND U */ /* STORED IN A/IW AFTER MA27B/BD. */ /* N - MUST BE SET TO THE ORDER OF THE MATRIX. NOT ALTERED */ /* BY MA27R/RD. */ /* A - MUST BE SET TO HOLD THE REAL VALUES CORRESPONDING */ /* TO THE FACTORS OF DINVERSE AND U. THIS MUST BE */ /* UNCHANGED SINCE THE PRECEDING CALL TO MA27B/BD. NOT ALTERED */ /* BY MA27R/RD. */ /* LA - LENGTH OF ARRAY A. NOT ALTERED BY MA27R/RD. */ /* IW - HOLDS THE INTEGER INDEXING */ /* INFORMATION FOR THE MATRIX FACTORS IN A. THIS MUST BE */ /* UNCHANGED SINCE THE PRECEDING CALL TO MA27B/BD. NOT ALTERED */ /* BY MA27R/RD. */ /* LIW - LENGTH OF ARRAY IW. NOT ALTERED BY MA27R/RD. */ /* W - USED */ /* AS WORKSPACE BY MA27R/RD TO HOLD THE COMPONENTS OF THE RIGHT */ /* HAND SIDES CORRESPONDING TO CURRENT BLOCK PIVOTAL ROWS. */ /* MAXFNT - INTEGER VARIABLE. MUST BE SET TO THE LARGEST NUMBER OF */ /* VARIABLES IN ANY BLOCK PIVOT ROW. THIS VALUE WAS GIVEN */ /* ON OUTPUT FROM MA27B/BD. NOT ALTERED BY MA27R/RD. */ /* RHS - ON INPUT, */ /* MUST BE SET TO HOLD THE RIGHT HAND SIDE MODIFIED BY THE */ /* FORWARD SUBSTITUTION OPERATIONS. ON OUTPUT, RHS WILL HOLD */ /* THE SOLUTION VECTOR. */ /* IW2 - ON ENTRY IW2(I) (I = 1,NBLK) */ /* MUST HOLD POINTERS TO THE */ /* BEGINNING OF EACH BLOCK PIVOT IN ARRAY IW, AS SET BY */ /* MA27Q/QD. */ /* NBLK - NUMBER OF BLOCK PIVOT ROWS. NOT ALTERED BY MA27R/RD. */ /* LATOP - IT IS THE POSITION IN */ /* A OF THE LAST ENTRY IN THE FACTORS. IT MUST BE UNCHANGED */ /* SINCE THE CALL TO MA27Q/QD. IT IS NOT ALTERED BY MA27R/RD. */ /* ICNTL is an INTEGER array of length 30, see MA27A/AD. */ /* ICNTL(IFRLVL+I) I=1,20 IS USED TO CONTROL WHETHER DIRECT OR INDIRECT */ /* ACCESS IS USED BY MA27C/CD. INDIRECT ACCESS IS EMPLOYED */ /* IN FORWARD AND BACK SUBSTITUTION RESPECTIVELY IF THE SIZE OF */ /* A BLOCK IS LESS THAN ICNTL(IFRLVL+MIN(10,NPIV)) AND */ /* ICNTL(IFRLVL+10+MIN(10,NPIV)) RESPECTIVELY, WHERE NPIV IS THE */ /* NUMBER OF PIVOTS IN THE BLOCK. */ /* .. Scalar Arguments .. */ /* .. */ /* .. Array Arguments .. */ /* .. */ /* .. Local Scalars .. */ /* .. */ /* .. Intrinsic Functions .. */ /* .. */ /* .. Executable Statements .. */ /* APOS. RUNNING POINTER TO CURRENT PIVOT POSITION IN ARRAY A. */ /* IPOS. RUNNING POINTER TO BEGINNING OF CURRENT BLOCK PIVOT ROW. */ /* Parameter adjustments */ --rhs; --a; --iw; --w; --iw2; --icntl; /* Function Body */ apos = *latop + 1; npiv = 0; iblk = *nblk + 1; /* RUN THROUGH BLOCK PIVOT ROWS IN THE REVERSE ORDER. */ i__1 = *n; for (loop = 1; loop <= i__1; ++loop) { if (npiv > 0) { goto L110; } --iblk; if (iblk < 1) { goto L190; } ipos = iw2[iblk]; /* ABS(LIELL) IS NUMBER OF VARIABLES (COLUMNS) IN BLOCK PIVOT ROW. */ liell = -iw[ipos]; /* NPIV IS NUMBER OF PIVOTS (ROWS) IN BLOCK PIVOT. */ npiv = 1; if (liell > 0) { goto L10; } liell = -liell; ++ipos; npiv = iw[ipos]; L10: jpos = ipos + npiv; j2 = ipos + liell; ilvl = min(10,npiv) + 10; if (liell < icntl[ilvl + 5]) { goto L110; } /* PERFORM OPERATIONS USING DIRECT ADDRESSING. */ j1 = ipos + 1; /* LOAD APPROPRIATE COMPONENTS OF RHS INTO W. */ ifr = 0; i__2 = j2; for (jj = j1; jj <= i__2; ++jj) { j = (i__3 = iw[jj], abs(i__3)); ++ifr; w[ifr] = rhs[j]; /* L20: */ } /* JPIV IS USED AS A FLAG SO THAT IPIV IS INCREMENTED CORRECTLY AFTER */ /* THE USE OF A 2 BY 2 PIVOT. */ jpiv = 1; /* PERFORM ELIMINATIONS. */ i__2 = npiv; for (iipiv = 1; iipiv <= i__2; ++iipiv) { --jpiv; if (jpiv == 1) { goto L90; } ipiv = npiv - iipiv + 1; if (ipiv == 1) { goto L30; } /* JUMP IF WE HAVE A 2 BY 2 PIVOT. */ if (iw[jpos - 1] < 0) { goto L60; } /* PERFORM BACK-SUBSTITUTION USING 1 BY 1 PIVOT. */ L30: jpiv = 1; apos -= liell + 1 - ipiv; ist = ipiv + 1; w1 = w[ipiv] * a[apos]; if (liell < ist) { goto L50; } jj1 = apos + 1; i__3 = liell; for (j = ist; j <= i__3; ++j) { w1 += a[jj1] * w[j]; ++jj1; /* L40: */ } L50: w[ipiv] = w1; --jpos; goto L90; /* PERFORM BACK-SUBSTITUTION OPERATIONS WITH 2 BY 2 PIVOT */ L60: jpiv = 2; apos2 = apos - (liell + 1 - ipiv); apos = apos2 - (liell + 2 - ipiv); ist = ipiv + 1; w1 = w[ipiv - 1] * a[apos] + w[ipiv] * a[apos + 1]; w2 = w[ipiv - 1] * a[apos + 1] + w[ipiv] * a[apos2]; if (liell < ist) { goto L80; } jj1 = apos + 2; jj2 = apos2 + 1; i__3 = liell; for (j = ist; j <= i__3; ++j) { w1 += w[j] * a[jj1]; w2 += w[j] * a[jj2]; ++jj1; ++jj2; /* L70: */ } L80: w[ipiv - 1] = w1; w[ipiv] = w2; jpos += -2; L90: ; } /* RELOAD WORKING VECTOR INTO SOLUTION VECTOR. */ ifr = 0; i__2 = j2; for (jj = j1; jj <= i__2; ++jj) { j = (i__3 = iw[jj], abs(i__3)); ++ifr; rhs[j] = w[ifr]; /* L100: */ } npiv = 0; goto L180; /* PERFORM OPERATIONS USING INDIRECT ADDRESSING. */ L110: if (npiv == 1) { goto L120; } /* JUMP IF WE HAVE A 2 BY 2 PIVOT. */ if (iw[jpos - 1] < 0) { goto L150; } /* PERFORM BACK-SUBSTITUTION USING 1 BY 1 PIVOT. */ L120: --npiv; apos -= j2 - jpos + 1; iirhs = iw[jpos]; w1 = rhs[iirhs] * a[apos]; j1 = jpos + 1; if (j1 > j2) { goto L140; } k = apos + 1; i__2 = j2; for (j = j1; j <= i__2; ++j) { irhs = (i__3 = iw[j], abs(i__3)); w1 += a[k] * rhs[irhs]; ++k; /* L130: */ } L140: rhs[iirhs] = w1; --jpos; goto L180; /* PERFORM OPERATIONS WITH 2 BY 2 PIVOT */ L150: npiv += -2; apos2 = apos - (j2 - jpos + 1); apos = apos2 - (j2 - jpos + 2); i1rhs = -iw[jpos - 1]; i2rhs = iw[jpos]; w1 = rhs[i1rhs] * a[apos] + rhs[i2rhs] * a[apos + 1]; w2 = rhs[i1rhs] * a[apos + 1] + rhs[i2rhs] * a[apos2]; j1 = jpos + 1; if (j1 > j2) { goto L170; } jj1 = apos + 2; jj2 = apos2 + 1; i__2 = j2; for (j = j1; j <= i__2; ++j) { irhs = (i__3 = iw[j], abs(i__3)); w1 += rhs[irhs] * a[jj1]; w2 += rhs[irhs] * a[jj2]; ++jj1; ++jj2; /* L160: */ } L170: rhs[i1rhs] = w1; rhs[i2rhs] = w2; jpos += -2; L180: ; } L190: return 0; } /* ma27r_ */