dsrc2c.c
来自「InsightToolkit-1.4.0(有大量的优化算法程序)」· C语言 代码 · 共 1,951 行 · 第 1/5 页
C
1,951 行
io___62.ciunit = itcom1_1.nout;
s_wsfe(&io___62);
do_fio(&c__1, (char *)&ier, (ftnlen)sizeof(integer));
e_wsfe();
}
goto L360;
L190:
if (itcom1_1.level <= 2) {
goto L210;
}
io___63.ciunit = itcom1_1.nout;
s_wsfe(&io___63);
e_wsfe();
L210:
if (iparm[11] != 0) {
goto L220;
}
timi1 = timer_(&dummy);
/* ... ITERATION SEQUENCE */
L220:
itmax1 = itcom1_1.itmax + 1;
i__1 = itmax1;
for (loop = 1; loop <= i__1; ++loop) {
itcom1_1.in = loop - 1;
if (itcom1_1.in % 2 == 1) {
goto L230;
}
/* ... CODE FOR THE EVEN ITERATIONS. */
/* U = U(IN) */
/* WKSP(IB1) = U(IN-1) */
itjsi_(&n, &ia[1], &ja[1], &a[1], &rhs[1], &u[1], &wksp[ib1], &wksp[
ib2], &icnt);
if (itcom2_1.halt) {
goto L270;
}
goto L240;
/* ... CODE FOR THE ODD ITERATIONS. */
/* U = U(IN-1) */
/* WKSP(IB1) = U(IN) */
L230:
itjsi_(&n, &ia[1], &ja[1], &a[1], &rhs[1], &wksp[ib1], &u[1], &wksp[
ib2], &icnt);
if (itcom2_1.halt) {
goto L270;
}
L240:
;
}
/* ... ITMAX HAS BEEN REACHED */
if (iparm[11] != 0) {
goto L250;
}
timi2 = timer_(&dummy);
time1 = (doublereal) (timi2 - timi1);
L250:
ier = 23;
if (itcom1_1.level >= 1) {
io___69.ciunit = itcom1_1.nout;
s_wsfe(&io___69);
do_fio(&c__1, (char *)&itcom1_1.itmax, (ftnlen)sizeof(integer));
e_wsfe();
}
if (iparm[3] == 0) {
rparm[1] = itcom3_1.stptst;
}
goto L300;
/* ... METHOD HAS CONVERGED */
L270:
if (iparm[11] != 0) {
goto L280;
}
timi2 = timer_(&dummy);
time1 = (doublereal) (timi2 - timi1);
L280:
if (itcom1_1.level >= 1) {
io___70.ciunit = itcom1_1.nout;
s_wsfe(&io___70);
do_fio(&c__1, (char *)&itcom1_1.in, (ftnlen)sizeof(integer));
e_wsfe();
}
/* ... PUT SOLUTION INTO U IF NOT ALREADY THERE. */
L300:
if (itcom1_1.in % 2 == 1) {
itpackdcopy_(&n, &wksp[ib1], &c__1, &u[1], &c__1);
}
/* ... UNSCALE THE MATRIX, SOLUTION, AND RHS VECTORS. */
unscal_(&n, &ia[1], &ja[1], &a[1], &rhs[1], &u[1], &wksp[1]);
/* ... UN-PERMUTE MATRIX,RHS, AND SOLUTION */
if (iparm[9] < 0) {
goto L330;
}
permat_(&n, &ia[1], &ja[1], &a[1], &iwksp[ib2], &iwksp[ib3], &
itcom1_1.isym, &itcom1_1.level, &itcom1_1.nout, &ierper);
if (ierper == 0) {
goto L320;
}
if (itcom1_1.level >= 0) {
io___72.ciunit = itcom1_1.nout;
s_wsfe(&io___72);
do_fio(&c__1, (char *)&ierper, (ftnlen)sizeof(integer));
e_wsfe();
}
if (ier == 0) {
ier = ierper;
}
goto L360;
L320:
pervec_(&n, &rhs[1], &iwksp[ib2]);
pervec_(&n, &u[1], &iwksp[ib2]);
/* ... OPTIONAL ERROR ANALYSIS */
L330:
idgts = iparm[12];
if (idgts < 0) {
goto L340;
}
if (iparm[2] <= 0) {
idgts = 0;
}
perror_(&n, &ia[1], &ja[1], &a[1], &rhs[1], &u[1], &wksp[1], &digit1, &
digit2, &idgts);
/* ... SET RETURN PARAMETERS IN IPARM AND RPARM */
L340:
if (iparm[11] != 0) {
goto L350;
}
timj2 = timer_(&dummy);
time2 = (doublereal) (timj2 - timj1);
L350:
if (iparm[3] != 0) {
goto L360;
}
iparm[1] = itcom1_1.in;
iparm[9] = nb;
rparm[2] = itcom3_1.cme;
rparm[3] = itcom3_1.sme;
rparm[9] = time1;
rparm[10] = time2;
rparm[11] = digit1;
rparm[12] = digit2;
L360:
*ierr = ier;
if (itcom1_1.level >= 3) {
echall_(&n, &ia[1], &ja[1], &a[1], &rhs[1], &iparm[1], &rparm[1], &
c__2);
}
return 0;
} /* jsi_ */
/* Subroutine */ int sor_(integer *nn, integer *ia, integer *ja, doublereal *
a, doublereal *rhs, doublereal *u, integer *iwksp, integer *nw,
doublereal *wksp, integer *iparm, doublereal *rparm, integer *ierr)
{
/* Format strings */
static char fmt_10[] = "(\0020\002///1x,\002BEGINNING OF ITPACK SOLUTION\
MODULE SOR\002)";
static char fmt_40[] = "(\0020\002,\002*** W A R N I N G ************\
\002/\0020\002,\002 IN ITPACK ROUTINE SOR\002/\002 \002,\002 RPARM(1) =\
\002,d10.3,\002 (ZETA)\002/\002 \002,\002 A VALUE THIS SMALL MAY HINDER C\
ONVERGENCE \002/\002 \002,\002 SINCE MACHINE PRECISION DRELPR =\002,d10.3/\
\002 \002,\002 ZETA RESET TO \002,d10.3)";
static char fmt_60[] = "(\0020\002,\002*** F A T A L E R R O R *****\
*******\002/\0020\002,\002 CALLED FROM ITPACK ROUTINE SOR \002/\002 \002\
,\002 INVALID MATRIX DIMENSION, N =\002,i8)";
static char fmt_80[] = "(\0020\002,\002*** F A T A L E R R O R *****\
*******\002/\0020\002,\002 CALLED FROM ITPACK ROUTINE SOR \002/\002 \002\
,\002 ERROR DETECTED IN SUBROUTINE SBELM \002/\002 \002,\002 WHICH RE\
MOVES ROWS AND COLUMNS OF SYSTEM \002/\002 \002,\002 WHEN DIAGONAL ENTRY \
TOO LARGE \002/\002 \002,\002 IER = \002,i5,5x,\002 RPARM(8) = \002,d10.\
3,\002 (TOL)\002)";
static char fmt_100[] = "(\0020\002,\002*** F A T A L E R R O R ****\
********\002/\0020\002,\002 CALLED FROM ITPACK ROUTINE SOR \002/\002 \002,\
\002 NOT ENOUGH WORKSPACE AT \002,i10/\002 \002,\002 SET IPARM(8) =\
\002,i10,\002 (NW)\002)";
static char fmt_120[] = "(\0020\002,\002*** F A T A L E R R O R ****\
********\002/\0020\002,\002 CALLED FROM ITPACK ROUTINE SOR \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 SOR \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 SOR \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,\002CME IS THE ESTIMATE OF THE LARGEST E\
IGENVALUE OF\002,\002 THE JACOBI MATRIX\002)";
static char fmt_210[] = "(1x,\002OMEGA IS THE RELAXATION FACTOR\002)";
static char fmt_260[] = "(\0020\002,\002*** W A R N I N G ***********\
*\002/\0020\002,\002 IN ITPACK ROUTINE SOR\002/\002 \002,\002 FAILURE \
TO CONVERGE IN\002,i5,\002 ITERATIONS\002)";
static char fmt_290[] = "(/1x,\002SOR 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 SOR \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 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 *);
extern doublereal timer_(real *);
static real dummy;
extern /* Subroutine */ int itsor_(integer *, integer *, integer *,
doublereal *, doublereal *, doublereal *, doublereal *);
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___75 = { 0, 0, 0, fmt_10, 0 };
static cilist io___81 = { 0, 0, 0, fmt_40, 0 };
static cilist io___86 = { 0, 0, 0, fmt_60, 0 };
static cilist io___88 = { 0, 0, 0, fmt_80, 0 };
static cilist io___92 = { 0, 0, 0, fmt_100, 0 };
static cilist io___95 = { 0, 0, 0, fmt_120, 0 };
static cilist io___96 = { 0, 0, 0, fmt_140, 0 };
static cilist io___97 = { 0, 0, 0, fmt_150, 0 };
static cilist io___98 = { 0, 0, 0, fmt_180, 0 };
static cilist io___99 = { 0, 0, 0, fmt_200, 0 };
static cilist io___100 = { 0, 0, 0, fmt_210, 0 };
static cilist io___105 = { 0, 0, 0, fmt_260, 0 };
static cilist io___106 = { 0, 0, 0, fmt_290, 0 };
static cilist io___108 = { 0, 0, 0, fmt_310, 0 };
/* ITPACK 2C MAIN SUBROUTINE SOR (SUCCESSIVE OVERRELATION) */
/* EACH OF THE MAIN SUBROUTINES: */
/* JCG, JSI, SOR, SSORCG, SSORSI, RSCG, RSSI */
/* CAN BE USED INDEPENDENTLY OF THE OTHERS */
/* ... FUNCTION: */
/* THIS SUBROUTINE, SOR, DRIVES THE SUCCESSIVE */
/* OVERRELAXATION 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. SOR NEEDS THIS */
/* TO BE IN LENGTH AT LEAST 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) */
/* ... SOR SUBPROGRAM REFERENCES: */
/* FROM ITPACK BISRCH, DFAULT, ECHALL, ECHOUT, IPSTR, ITERM, */
/* TIMER, ITSOR, IVFILL, PERMAT, PERROR, */
/* PERVEC, PFSOR1, PMULT, PRBNDX, PSTOP, QSORT, */
/* SBELM, SCAL, DCOPY, DDOT, TAU, UNSCAL, VFILL, */
/* VOUT, WEVMW */
/* SYSTEM DABS, DLOG10, DBLE(AMAX0), DMAX1, DBLE(FLOAT), */
/* 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 */
⌨️ 快捷键说明
复制代码Ctrl + C
搜索代码Ctrl + F
全屏模式F11
增大字号Ctrl + =
减小字号Ctrl + -
显示快捷键?