dsrc2c.c
来自「InsightToolkit-1.4.0(有大量的优化算法程序)」· C语言 代码 · 共 1,951 行 · 第 1/5 页
C
1,951 行
********\002/\0020\002,\002 CALLED FROM ITPACK ROUTINE JSI \002/\002 \002,\
\002 ERROR DETECTED IN SUBROUTINE PRBNDX\002/\002 \002,\002 WHICH COM\
PUTES THE RED-BLACK INDEXING\002/\002 \002,\002 IER = \002,i5,\002 IPARM(\
9) = \002,i5,\002 (NB)\002)";
static char fmt_140[] = "(/10x,\002ORDER OF BLACK SUBSYSTEM = \002,i5\
,\002 (NB)\002)";
static char fmt_150[] = "(\0020\002,\002*** F A T A L E R R O R ****\
********\002/\0020\002,\002 CALLED FROM ITPACK ROUTINE JSI \002/\002 \002,\
\002 ERROR DETECTED IN SUBROUTINE PERMAT\002/\002 \002,\002 WHICH DOE\
S THE RED-BLACK PERMUTATION\002/\002 \002,\002 IER = \002,i5)";
static char fmt_180[] = "(\0020\002,\002*** F A T A L E R R O R ****\
********\002/\0020\002,\002 CALLED FROM ITPACK ROUTINE JSI \002/\002 \002,\
\002 ERROR DETECTED IN SUBROUTINE SCAL \002/\002 \002,\002 WHICH SCA\
LES THE SYSTEM \002/\002 \002,\002 IER = \002,i5)";
static char fmt_200[] = "(///1x,\002IN THE FOLLOWING, RHO AND GAMMA AR\
E\002,\002 ACCELERATION PARAMETERS\002)";
static char fmt_260[] = "(\0020\002,\002*** W A R N I N G ***********\
*\002/\0020\002,\002 IN ITPACK ROUTINE JSI\002/\002 \002,\002 FAILURE \
TO CONVERGE IN\002,i5,\002 ITERATIONS\002)";
static char fmt_290[] = "(/1x,\002JSI HAS CONVERGED IN \002,i5,\002 ITE\
RATIONS\002)";
static char fmt_310[] = "(\0020\002,\002*** F A T A L E R R O R ****\
********\002/\0020\002,\002 CALLED FROM ITPACK ROUTINE JSI \002/\002 \002,\
\002 ERROR DETECTED IN SUBROUTINE PERMAT\002/\002 \002,\002 WHICH UND\
OES THE RED-BLACK PERMUTATION \002/\002 \002,\002 IER = \002,i5)";
/* System generated locals */
integer i__1;
/* Builtin functions */
integer s_wsfe(cilist *), e_wsfe(void), do_fio(integer *, char *, ftnlen);
/* Local variables */
static integer n, n3, nb, ib1, ib2, ib3, ier;
static doublereal tol;
extern /* Subroutine */ int scal_(integer *, integer *, integer *,
doublereal *, doublereal *, doublereal *, doublereal *, integer *,
integer *, integer *);
static integer icnt;
static doublereal temp;
static integer loop;
static doublereal time1, time2;
static real timi1, timj1, timi2, timj2;
extern /* Subroutine */ int sbelm_(integer *, integer *, integer *,
doublereal *, doublereal *, integer *, doublereal *, doublereal *,
integer *, integer *, integer *, integer *);
static integer idgts;
extern /* Subroutine */ int vfill_(integer *, doublereal *, doublereal *),
itpackdcopy_(integer *, doublereal *, integer *, doublereal *, integer
*);
extern doublereal timer_(real *);
extern /* Subroutine */ int itjsi_(integer *, integer *, integer *,
doublereal *, doublereal *, doublereal *, doublereal *,
doublereal *, integer *);
static real dummy;
static doublereal digit1, digit2;
static integer itmax1;
extern /* Subroutine */ int echall_(integer *, integer *, integer *,
doublereal *, doublereal *, integer *, doublereal *, integer *),
pervec_(integer *, doublereal *, integer *), ivfill_(integer *,
integer *, integer *);
static integer ierper;
extern /* Subroutine */ int echout_(integer *, doublereal *, integer *),
permat_(integer *, integer *, integer *, doublereal *, integer *,
integer *, integer *, integer *, integer *, integer *), unscal_(
integer *, integer *, integer *, doublereal *, doublereal *,
doublereal *, doublereal *), prbndx_(integer *, integer *,
integer *, integer *, integer *, integer *, integer *, integer *,
integer *), perror_(integer *, integer *, integer *, doublereal *,
doublereal *, doublereal *, doublereal *, doublereal *,
doublereal *, integer *);
/* Fortran I/O blocks */
static cilist io___39 = { 0, 0, 0, fmt_10, 0 };
static cilist io___45 = { 0, 0, 0, fmt_40, 0 };
static cilist io___50 = { 0, 0, 0, fmt_60, 0 };
static cilist io___52 = { 0, 0, 0, fmt_80, 0 };
static cilist io___56 = { 0, 0, 0, fmt_100, 0 };
static cilist io___59 = { 0, 0, 0, fmt_120, 0 };
static cilist io___60 = { 0, 0, 0, fmt_140, 0 };
static cilist io___61 = { 0, 0, 0, fmt_150, 0 };
static cilist io___62 = { 0, 0, 0, fmt_180, 0 };
static cilist io___63 = { 0, 0, 0, fmt_200, 0 };
static cilist io___69 = { 0, 0, 0, fmt_260, 0 };
static cilist io___70 = { 0, 0, 0, fmt_290, 0 };
static cilist io___72 = { 0, 0, 0, fmt_310, 0 };
/* ITPACK 2C MAIN SUBROUTINE JSI (JACOBI SEMI-ITERATIVE) */
/* EACH OF THE MAIN SUBROUTINES: */
/* JCG, JSI, SOR, SSORCG, SSORSI, RSCG, RSSI */
/* CAN BE USED INDEPENDENTLY OF THE OTHERS */
/* ... FUNCTION: */
/* THIS SUBROUTINE, JSI, DRIVES THE JACOBI SEMI- */
/* ITERATION ALGORITHM. */
/* ... PARAMETER LIST: */
/* N INPUT INTEGER. DIMENSION OF THE MATRIX. (= NN) */
/* IA,JA INPUT INTEGER VECTORS. THE TWO INTEGER ARRAYS OF */
/* THE SPARSE MATRIX REPRESENTATION. */
/* A INPUT D.P. VECTOR. THE D.P. ARRAY OF THE SPARSE */
/* MATRIX REPRESENTATION. */
/* RHS INPUT D.P. VECTOR. CONTAINS THE RIGHT HAND SIDE */
/* OF THE MATRIX PROBLEM. */
/* U INPUT/OUTPUT D.P. VECTOR. ON INPUT, U CONTAINS THE */
/* INITIAL GUESS TO THE SOLUTION. ON OUTPUT, IT CONTAINS */
/* THE LATEST ESTIMATE TO THE SOLUTION. */
/* IWKSP INTEGER VECTOR WORKSPACE OF LENGTH 3*N */
/* NW INPUT INTEGER. LENGTH OF AVAILABLE WKSP. ON OUTPUT, */
/* IPARM(8) IS AMOUNT USED. */
/* WKSP D.P. VECTOR USED FOR WORKING SPACE. JACOBI SI */
/* NEEDS THIS TO BE IN LENGTH AT LEAST */
/* 2*N */
/* IPARM INTEGER VECTOR OF LENGTH 12. ALLOWS USER TO SPECIFY */
/* SOME INTEGER PARAMETERS WHICH AFFECT THE METHOD. */
/* RPARM D.P. VECTOR OF LENGTH 12. ALLOWS USER TO SPECIFY SOME */
/* D.P. PARAMETERS WHICH AFFECT THE METHOD. */
/* IER OUTPUT INTEGER. ERROR FLAG. (= IERR) */
/* ... JSI SUBPROGRAM REFERENCES: */
/* FROM ITPACK BISRCH, CHEBY, CHGSI, CHGSME, DFAULT, ECHALL, */
/* ECHOUT, ITERM, TIMER, ITJSI, IVFILL, PAR */
/* PERMAT, PERROR, PERVEC, PJAC, PMULT, PRBNDX, */
/* PSTOP, PVTBV, QSORT, DAXPY, SBELM, SCAL, */
/* DCOPY, DDOT, SUM3, TSTCHG, UNSCAL, VEVMW, */
/* VFILL, VOUT, WEVMW */
/* SYSTEM DABS, DLOG10, DBLE(AMAX0), DMAX1, DBLE(FLOAT), */
/* MOD,DSQRT */
/* VERSION: ITPACK 2C (MARCH 1982) */
/* CODE WRITTEN BY: DAVID KINCAID, ROGER GRIMES, JOHN RESPESS */
/* CENTER FOR NUMERICAL ANALYSIS */
/* UNIVERSITY OF TEXAS */
/* AUSTIN, TX 78712 */
/* (512) 471-1242 */
/* FOR ADDITIONAL DETAILS ON THE */
/* (A) SUBROUTINE SEE TOMS ARTICLE 1982 */
/* (B) ALGORITHM SEE CNA REPORT 150 */
/* BASED ON THEORY BY: DAVID YOUNG, DAVID KINCAID, LOU HAGEMAN */
/* REFERENCE THE BOOK: APPLIED ITERATIVE METHODS */
/* L. HAGEMAN, D. YOUNG */
/* ACADEMIC PRESS, 1981 */
/* ************************************************** */
/* * IMPORTANT NOTE * */
/* * * */
/* * WHEN INSTALLING ITPACK ROUTINES ON A * */
/* * DIFFERENT COMPUTER, RESET SOME OF THE VALUES * */
/* * IN SUBROUTNE DFAULT. MOST IMPORTANT ARE * */
/* * * */
/* * DRELPR MACHINE RELATIVE PRECISION * */
/* * RPARM(1) STOPPING CRITERION * */
/* * * */
/* * ALSO CHANGE SYSTEM-DEPENDENT ROUTINE * */
/* * SECOND USED IN TIMER * */
/* * * */
/* ************************************************** */
/* SPECIFICATIONS FOR ARGUMENTS */
/* SPECIFICATIONS FOR LOCAL VARIABLES */
/* *** BEGIN: ITPACK COMMON */
/* *** END : ITPACK COMMON */
/* ... VARIABLES IN COMMON BLOCK - ITCOM1 */
/* IN - ITERATION NUMBER */
/* IS - ITERATION NUMBER WHEN PARAMETERS LAST CHANGED */
/* ISYM - SYMMETRIC/NONSYMMETRIC STORAGE FORMAT SWITCH */
/* ITMAX - MAXIMUM NUMBER OF ITERATIONS ALLOWED */
/* LEVEL - LEVEL OF OUTPUT CONTROL SWITCH */
/* NOUT - OUTPUT UNIT NUMBER */
/* ... VARIABLES IN COMMON BLOCK - ITCOM2 */
/* ADAPT - FULLY ADAPTIVE PROCEDURE SWITCH */
/* BETADT - SWITCH FOR ADAPTIVE DETERMINATION OF BETA */
/* CASEII - ADAPTIVE PROCEDURE CASE SWITCH */
/* HALT - STOPPING TEST SWITCH */
/* PARTAD - PARTIALLY ADAPTIVE PROCEDURE SWITCH */
/* ... VARIABLES IN COMMON BLOCK - ITCOM3 */
/* BDELNM - TWO NORM OF B TIMES DELTA-SUPER-N */
/* BETAB - ESTIMATE FOR THE SPECTRAL RADIUS OF LU MATRIX */
/* CME - ESTIMATE OF LARGEST EIGENVALUE */
/* DELNNM - INNER PRODUCT OF PSEUDO-RESIDUAL AT ITERATION N */
/* DELSNM - INNER PRODUCT OF PSEUDO-RESIDUAL AT ITERATION S */
/* FF - ADAPTIVE PROCEDURE DAMPING FACTOR */
/* GAMMA - ACCELERATION PARAMETER */
/* OMEGA - OVERRELAXATION PARAMETER FOR SOR AND SSOR */
/* QA - PSEUDO-RESIDUAL RATIO */
/* QT - VIRTUAL SPECTRAL RADIUS */
/* RHO - ACCELERATION PARAMETER */
/* RRR - ADAPTIVE PARAMETER */
/* SIGE - PARAMETER SIGMA-SUB-E */
/* SME - ESTIMATE OF SMALLEST EIGENVALUE */
/* SPECR - SPECTRAL RADIUS ESTIMATE FOR SSOR */
/* DRELPR - MACHINE RELATIVE PRECISION */
/* STPTST - STOPPING PARAMETER */
/* UDNM - TWO NORM OF U */
/* ZETA - STOPPING CRITERION */
/* ... INITIALIZE COMMON BLOCKS */
/* Parameter adjustments */
--u;
--rhs;
--ia;
--ja;
--a;
--iwksp;
--wksp;
--iparm;
--rparm;
/* Function Body */
itcom1_1.level = iparm[2];
itcom1_1.nout = iparm[4];
if (itcom1_1.level >= 1) {
io___39.ciunit = itcom1_1.nout;
s_wsfe(&io___39);
e_wsfe();
}
ier = 0;
if (iparm[1] <= 0) {
return 0;
}
n = *nn;
if (iparm[11] == 0) {
timj1 = timer_(&dummy);
}
if (itcom1_1.level >= 3) {
goto L20;
}
echout_(&iparm[1], &rparm[1], &c__2);
goto L30;
L20:
echall_(&n, &ia[1], &ja[1], &a[1], &rhs[1], &iparm[1], &rparm[1], &c__1);
L30:
temp = itcom3_1.drelpr * 500.;
if (itcom3_1.zeta >= temp) {
goto L50;
}
if (itcom1_1.level >= 1) {
io___45.ciunit = itcom1_1.nout;
s_wsfe(&io___45);
do_fio(&c__1, (char *)&itcom3_1.zeta, (ftnlen)sizeof(doublereal));
do_fio(&c__1, (char *)&itcom3_1.drelpr, (ftnlen)sizeof(doublereal));
do_fio(&c__1, (char *)&temp, (ftnlen)sizeof(doublereal));
e_wsfe();
}
itcom3_1.zeta = temp;
L50:
time1 = rparm[9];
time2 = rparm[10];
digit1 = rparm[11];
digit2 = rparm[12];
/* ... VERIFY N */
if (n > 0) {
goto L70;
}
ier = 21;
if (itcom1_1.level >= 0) {
io___50.ciunit = itcom1_1.nout;
s_wsfe(&io___50);
do_fio(&c__1, (char *)&n, (ftnlen)sizeof(integer));
e_wsfe();
}
goto L360;
L70:
/* ... REMOVE ROWS AND COLUMNS IF REQUESTED */
if (iparm[10] == 0) {
goto L90;
}
tol = rparm[8];
ivfill_(&n, &iwksp[1], &c__0);
vfill_(&n, &wksp[1], &c_b21);
sbelm_(&n, &ia[1], &ja[1], &a[1], &rhs[1], &iwksp[1], &wksp[1], &tol, &
itcom1_1.isym, &itcom1_1.level, &itcom1_1.nout, &ier);
if (ier == 0) {
goto L90;
}
if (itcom1_1.level >= 0) {
io___52.ciunit = itcom1_1.nout;
s_wsfe(&io___52);
do_fio(&c__1, (char *)&ier, (ftnlen)sizeof(integer));
do_fio(&c__1, (char *)&tol, (ftnlen)sizeof(doublereal));
e_wsfe();
}
goto L360;
/* ... INITIALIZE WKSP BASE ADDRESSES. */
L90:
ib1 = 1;
ib2 = ib1 + n;
ib3 = ib2 + n;
iparm[8] = n << 1;
if (*nw >= iparm[8]) {
goto L110;
}
ier = 22;
if (itcom1_1.level >= 0) {
io___56.ciunit = itcom1_1.nout;
s_wsfe(&io___56);
do_fio(&c__1, (char *)&(*nw), (ftnlen)sizeof(integer));
do_fio(&c__1, (char *)&iparm[8], (ftnlen)sizeof(integer));
e_wsfe();
}
goto L360;
/* ... PERMUTE TO RED-BLACK SYSTEM IF REQUESTED */
L110:
nb = iparm[9];
if (nb < 0) {
goto L170;
}
n3 = n * 3;
ivfill_(&n3, &iwksp[1], &c__0);
prbndx_(&n, &nb, &ia[1], &ja[1], &iwksp[1], &iwksp[ib2], &itcom1_1.level,
&itcom1_1.nout, &ier);
if (ier == 0) {
goto L130;
}
if (itcom1_1.level >= 0) {
io___59.ciunit = itcom1_1.nout;
s_wsfe(&io___59);
do_fio(&c__1, (char *)&ier, (ftnlen)sizeof(integer));
do_fio(&c__1, (char *)&nb, (ftnlen)sizeof(integer));
e_wsfe();
}
goto L360;
/* ... PERMUTE MATRIX AND RHS */
L130:
if (itcom1_1.level >= 2) {
io___60.ciunit = itcom1_1.nout;
s_wsfe(&io___60);
do_fio(&c__1, (char *)&nb, (ftnlen)sizeof(integer));
e_wsfe();
}
permat_(&n, &ia[1], &ja[1], &a[1], &iwksp[1], &iwksp[ib3], &itcom1_1.isym,
&itcom1_1.level, &itcom1_1.nout, &ier);
if (ier == 0) {
goto L160;
}
if (itcom1_1.level >= 0) {
io___61.ciunit = itcom1_1.nout;
s_wsfe(&io___61);
do_fio(&c__1, (char *)&ier, (ftnlen)sizeof(integer));
e_wsfe();
}
goto L360;
L160:
pervec_(&n, &rhs[1], &iwksp[1]);
pervec_(&n, &u[1], &iwksp[1]);
/* ... SCALE LINEAR SYSTEM, U, AND RHS BY THE SQUARE ROOT OF THE */
/* ... DIAGONAL ELEMENTS. */
L170:
vfill_(&iparm[8], &wksp[1], &c_b21);
scal_(&n, &ia[1], &ja[1], &a[1], &rhs[1], &u[1], &wksp[1], &
itcom1_1.level, &itcom1_1.nout, &ier);
if (ier == 0) {
goto L190;
}
if (itcom1_1.level >= 0) {
⌨️ 快捷键说明
复制代码Ctrl + C
搜索代码Ctrl + F
全屏模式F11
增大字号Ctrl + =
减小字号Ctrl + -
显示快捷键?