📄 dsrc2c.c
字号:
/* ************************************************** */
/* */
/* SPECIFICATIONS FOR ARGUMENTS */
/* */
/* SPECIFICATIONS FOR LOCAL VARIABLES */
/* */
/* ... 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 */
itcom1_1.level = iparm[1];
itcom1_1.nout = iparm[3];
ier = 0;
if (iparm[0] <= 0)
return 0;
if (iparm[10] == 0)
timj1 = dsrc_timer_((real*)0);
if (itcom1_1.level < 3)
echout_(iparm, rparm, &c__2);
else
echall_(n, ia, ja, a, rhs, iparm, rparm, &c__1);
temp = itcom3_1.drelpr * 500.;
if (itcom3_1.zeta < temp)
itcom3_1.zeta = temp;
time1 = rparm[8];
time2 = rparm[9];
digit1 = rparm[10];
digit2 = rparm[11];
/* ... VERIFY N */
if (*n <= 0) {
ier = 21;
goto L360;
}
/* ... REMOVE ROWS AND COLUMNS IF REQUESTED */
if (iparm[9] != 0) {
tol = rparm[7];
ivfill_(n, iwksp, &c__0);
vfill_(n, wksp, &c_b21);
sbelm_(n, ia, ja, a, rhs, iwksp, wksp, &tol, &itcom1_1.isym, &itcom1_1.level, &itcom1_1.nout, &ier);
if (ier != 0)
goto L360;
}
/* ... INITIALIZE WKSP BASE ADDRESSES. */
ib1 = 0;
ib2 = ib1 + *n;
ib3 = ib2 + *n;
iparm[7] = *n << 1;
if (*nw < iparm[7]) {
ier = 22;
goto L360;
}
/* ... PERMUTE TO RED-BLACK SYSTEM IF REQUESTED */
nb = iparm[8];
if (nb < 0)
goto L170;
n3 = *n * 3;
ivfill_(&n3, iwksp, &c__0);
prbndx_(n, &nb, ia, ja, iwksp, &iwksp[ib2], &itcom1_1.level, &itcom1_1.nout, &ier);
if (ier != 0)
goto L360;
/* ... PERMUTE MATRIX AND RHS */
permat_(n, ia, ja, a, iwksp, &iwksp[ib3], &itcom1_1.isym, &itcom1_1.level, &itcom1_1.nout, &ier);
if (ier != 0)
goto L360;
pervec_(n, rhs, iwksp);
pervec_(n, u, iwksp);
/* ... SCALE LINEAR SYSTEM, U, AND RHS BY THE SQUARE ROOT OF THE DIAGONAL ELEMENTS. */
L170:
vfill_(&iparm[7], wksp, &c_b21);
scal_(n, ia, ja, a, rhs, u, wksp, &itcom1_1.level, &itcom1_1.nout, &ier);
if (ier != 0)
goto L360;
if (iparm[10] == 0)
timi1 = dsrc_timer_((real*)0);
/* ... ITERATION SEQUENCE */
itmax1 = itcom1_1.itmax + 1;
for (loop = 1; loop <= itmax1; ++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, ja, a, rhs, u, &wksp[ib1], &wksp[ib2], &icnt);
if (itcom2_1.halt)
goto L270;
continue;
/* ... CODE FOR THE ODD ITERATIONS. */
/* U = U(IN-1) */
/* WKSP(IB1) = U(IN) */
L230:
itjsi_(n, ia, ja, a, rhs, &wksp[ib1], u, &wksp[ib2], &icnt);
if (itcom2_1.halt)
goto L270;
}
/* ... ITMAX HAS BEEN REACHED */
if (iparm[10] == 0) {
timi2 = dsrc_timer_((real*)0);
time1 = (doublereal) (timi2 - timi1);
}
ier = 23;
if (iparm[2] == 0)
rparm[0] = itcom3_1.stptst;
goto L300;
/* ... METHOD HAS CONVERGED */
L270:
if (iparm[10] == 0) {
timi2 = dsrc_timer_((real*)0);
time1 = (doublereal) (timi2 - timi1);
}
/* ... PUT SOLUTION INTO U IF NOT ALREADY THERE. */
L300:
if (itcom1_1.in % 2 == 1)
itpackdcopy_(n, &wksp[ib1], &c__1, u, &c__1);
/* ... UNSCALE THE MATRIX, SOLUTION, AND RHS VECTORS. */
unscal_(n, ia, ja, a, rhs, u, wksp);
/* ... UN-PERMUTE MATRIX,RHS, AND SOLUTION */
if (iparm[8] < 0)
goto L330;
permat_(n, ia, ja, a, &iwksp[ib2], &iwksp[ib3], &itcom1_1.isym, &itcom1_1.level, &itcom1_1.nout, &ierper);
if (ierper != 0) {
if (ier == 0)
ier = ierper;
goto L360;
}
pervec_(n, rhs, &iwksp[ib2]);
pervec_(n, u, &iwksp[ib2]);
/* ... OPTIONAL ERROR ANALYSIS */
L330:
idgts = iparm[11];
if (idgts >= 0) {
if (iparm[1] <= 0)
idgts = 0;
perror_(n, ia, ja, a, rhs, u, wksp, &digit1, &digit2, &idgts);
}
/* ... SET RETURN PARAMETERS IN IPARM AND RPARM */
if (iparm[10] == 0) {
timj2 = dsrc_timer_((real*)0);
time2 = (doublereal) (timj2 - timj1);
}
if (iparm[2] == 0) {
iparm[0] = itcom1_1.in;
iparm[8] = nb;
rparm[1] = itcom3_1.cme;
rparm[2] = itcom3_1.sme;
rparm[8] = time1;
rparm[9] = time2;
rparm[10] = digit1;
rparm[11] = digit2;
}
L360:
*ierr = ier;
if (itcom1_1.level >= 3)
echall_(n, ia, ja, a, rhs, iparm, rparm, &c__2);
return 0;
} /* jsi_ */
/* Subroutine */
int sor_(integer *n, integer *ia, integer *ja, doublereal *a, doublereal *rhs, doublereal *u,
integer *iwksp, integer *nw, doublereal *wksp, integer *iparm, doublereal *rparm, integer* ierr)
{
/* Local variables */
static integer n3, nb, ib1, ib2, ib3, ier;
static doublereal tol;
static doublereal temp;
static integer loop;
static doublereal time1, time2;
static real timi1, timj1, timi2, timj2;
static integer idgts;
static doublereal digit1, digit2;
static integer itmax1;
static integer ierper;
/* 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 */
/* THIS SUBROUTINE, SOR, DRIVES THE SUCCESSIVE */
/* OVERRELAXATION ALGORITHM. */
/* */
/* ... PARAMETER LIST: */
/* */
/* N INPUT INTEGER. DIMENSION OF THE MATRIX. */
/* 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) */
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -