📄 dsrc2c.c
字号:
/* */
/* 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];
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__3);
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 = 31;
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;
if (*nw < iparm[7]) {
ier = 32;
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;
/* ... CODE FOR ONE ITERATION. */
/* U = U(IN) */
itsor_(n, ia, ja, a, rhs, u, &wksp[ib1]);
if (itcom2_1.halt)
goto L270;
}
/* ... ITMAX HAS BEEN REACHED */
if (iparm[10] == 0) {
timi2 = dsrc_timer_((real*)0);
time1 = (doublereal) (timi2 - timi1);
}
ier = 33;
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);
}
/* ... UNSCALE THE MATRIX, SOLUTION, AND RHS VECTORS. */
L300:
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[4] = itcom3_1.omega;
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;
} /* sor_ */
/* Subroutine */
int ssorcg_(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, ib4, ib5, ib6, ib7, 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 doublereal betnew;
static integer ierper;
/* ITPACK 2C MAIN SUBROUTINE SSORCG (SYMMETRIC SUCCESSIVE OVER- */
/* RELAXATION CONJUGATE GRADIENT) */
/* EACH OF THE MAIN SUBROUTINES: */
/* JCG, JSI, SOR, SSORCG, SSORSI, RSCG, RSSI */
/* CAN BE USED INDEPENDENTLY OF THE OTHERS */
/* THIS SUBROUTINE, SSORCG, DRIVES THE SYMMETRIC SOR-CG */
/* 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. SSOR-CG */
/* NEEDS TO BE IN LENGTH AT LEAST */
/* 6*N + 2*ITMAX, IF IPARM(5)=0 (SYMMETRIC STORAGE) */
/* 6*N + 4*ITMAX, IF IPARM(5)=1 (NONSYMMETRIC STORAGE) */
/* IPARM INTEGER VECTOR OF LENGTH 12. ALLOWS USER TO SPECIFY */
/* SOME INTEGER PARAMETERS WHICH AFFECT THE METHOD. IF */
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -