📄 dsrc2c.c
字号:
/* 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) */
/* */
/* ... SSORCG SUBPROGRAM REFERENCES: */
/* */
/* FROM ITPACK BISRCH, CHGCON, DETERM, DFAULT, ECHALL, */
/* ECHOUT, EIGVNS, EIGVSS, EQRT1S, ITERM, TIMER, */
/* ITSRCG, IVFILL, OMEG, OMGCHG, OMGSTR, */
/* PARCON, PBETA, PBSOR, PERMAT, PERROR, */
/* PERVEC, PFSOR, PJAC, PMULT, PRBNDX, PSTOP, PVT */
/* QSORT, SBELM, SCAL, DCOPY, DDOT, SUM3, */
/* UNSCAL, VEVMW, VEVPW, VFILL, VOUT, WEVMW, */
/* ZBRENT */
/* SYSTEM DABS, DLOG, DLOG10, DBLE(AMAX0), DMAX1, AMIN1, */
/* 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 */
/* */
/* ... 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];
if (iparm[8] >= 0)
iparm[5] = 2;
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__4);
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 = 41;
goto L390;
}
/* ... 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 L390;
}
/* ... INITIALIZE WKSP BASE ADDRESSES. */
ib1 = 0;
ib2 = ib1 + *n;
ib3 = ib2 + *n;
ib4 = ib3 + *n;
ib5 = ib4 + *n;
ib6 = ib5 + *n;
ib7 = ib6 + *n;
iparm[7] = *n * 6 + (itcom1_1.itmax << 1);
if (itcom1_1.isym != 0)
iparm[7] += itcom1_1.itmax << 1;
if (*nw < iparm[7]) {
ier = 42;
goto L390;
}
/* ... 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 L390;
/* ... 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 L390;
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 L390;
if (iparm[10] == 0)
timi1 = dsrc_timer_((real*)0);
/* ... SPECIAL PROCEDURE FOR FULLY ADAPTIVE CASE. */
if (! itcom2_1.adapt)
goto L250;
if (itcom2_1.betadt) {
vfill_(n, &wksp[ib1], &c_b286);
betnew = pbeta_(n, ia, ja, a, &wksp[ib1], &wksp[ib2], &wksp[ib3]) / (doublereal)(*n);
itcom3_1.betab = max(max(itcom3_1.betab,.25),betnew);
}
omeg_(&c_b21, &c__1);
itcom1_1.is = 0;
/* ... INITIALIZE FORWARD PSEUDO-RESIDUAL */
L250:
itpackdcopy_(n, rhs, &c__1, &wksp[ib1], &c__1);
itpackdcopy_(n, u, &c__1, &wksp[ib2], &c__1);
pfsor_(n, ia, ja, a, &wksp[ib2], &wksp[ib1]);
vevmw_(n, &wksp[ib2], u);
/* ... 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 L260;
/* ... CODE FOR THE EVEN ITERATIONS. */
/* U = U(IN) WKSP(IB2) = C(IN) */
/* WKSP(IB1) = U(IN-1) WKSP(IB3) = C(IN-1) */
itsrcg_(n, ia, ja, a, rhs, u, &wksp[ib1], &wksp[ib2], &wksp[ib3], &wksp[ib4], &wksp[ib5], &wksp[ib6], &wksp[ib7]);
if (itcom2_1.halt)
goto L300;
continue;
/* ... CODE FOR THE ODD ITERATIONS. */
/* U = U(IN-1) WKSP(IB2) = C(IN-1) */
/* WKSP(IB1) = U(IN) WKSP(IB3) =C(IN) */
L260:
itsrcg_(n, ia, ja, a, rhs, &wksp[ib1], u, &wksp[ib3], &wksp[ib2], &wksp[ib4], &wksp[ib5], &wksp[ib6], &wksp[ib7]);
if (itcom2_1.halt)
goto L300;
}
/* ... ITMAX HAS BEEN REACHED */
if (iparm[10] == 0) {
timi2 = dsrc_timer_((real*)0);
time1 = (doublereal) (timi2 - timi1);
}
ier = 43;
if (iparm[2] == 0)
rparm[0] = itcom3_1.stptst;
goto L330;
/* ... METHOD HAS CONVERGED */
L300:
if (iparm[10] == 0) {
timi2 = dsrc_timer_((real*)0);
time1 = (doublereal) (timi2 - timi1);
}
/* ... PUT SOLUTION INTO U IF NOT ALREADY THERE. */
L330:
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 L360;
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 L390;
}
pervec_(n, rhs, &iwksp[ib2]);
pervec_(n, u, &iwksp[ib2]);
/* ... OPTIONAL ERROR ANALYSIS */
L360:
idgts = iparm[11];
if (idgts >= 0) {
if (iparm[1] <= 0)
idgts = 0;
perror_(n, ia, ja, a, rhs, u, wksp, &digit1, &digit2, &idgts);
}
/* ... SET⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -