dsrc2c.c

来自「InsightToolkit-1.4.0(有大量的优化算法程序)」· C语言 代码 · 共 1,951 行 · 第 1/5 页

/*     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___75.ciunit = itcom1_1.nout;
  s_wsfe(&io___75);
  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__3);
    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___81.ciunit = itcom1_1.nout;
  s_wsfe(&io___81);
  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 = 31;
    if (itcom1_1.level >= 0) {
  io___86.ciunit = itcom1_1.nout;
  s_wsfe(&io___86);
  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___88.ciunit = itcom1_1.nout;
  s_wsfe(&io___88);
  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;
    if (*nw >= iparm[8]) {
  goto L110;
    }
    ier = 32;
    if (itcom1_1.level >= 0) {
  io___92.ciunit = itcom1_1.nout;
  s_wsfe(&io___92);
  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___95.ciunit = itcom1_1.nout;
  s_wsfe(&io___95);
  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___96.ciunit = itcom1_1.nout;
  s_wsfe(&io___96);
  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___97.ciunit = itcom1_1.nout;
  s_wsfe(&io___97);
  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) {
  io___98.ciunit = itcom1_1.nout;
  s_wsfe(&io___98);
  do_fio(&c__1, (char *)&ier, (ftnlen)sizeof(integer));
  e_wsfe();
    }
    goto L360;
L190:
    if (itcom1_1.level <= 2) {
  goto L220;
    }
    if (itcom2_1.adapt) {
  io___99.ciunit = itcom1_1.nout;
  s_wsfe(&io___99);
  e_wsfe();
    }
    io___100.ciunit = itcom1_1.nout;
    s_wsfe(&io___100);
    e_wsfe();
L220:
    if (iparm[11] != 0) {
  goto L230;
    }
    timi1 = timer_(&dummy);

/* ... ITERATION SEQUENCE */

L230:
    itmax1 = itcom1_1.itmax + 1;
    i__1 = itmax1;
    for (loop = 1; loop <= i__1; ++loop) {
  itcom1_1.in = loop - 1;

/* ... CODE FOR ONE ITERATION. */

/*     U           = U(IN) */

  itsor_(&n, &ia[1], &ja[1], &a[1], &rhs[1], &u[1], &wksp[ib1]);

  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:
    if (itcom1_1.level >= 1) {
  io___105.ciunit = itcom1_1.nout;
  s_wsfe(&io___105);
  do_fio(&c__1, (char *)&itcom1_1.itmax, (ftnlen)sizeof(integer));
  e_wsfe();
    }
    ier = 33;
    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___106.ciunit = itcom1_1.nout;
  s_wsfe(&io___106);
  do_fio(&c__1, (char *)&itcom1_1.in, (ftnlen)sizeof(integer));
  e_wsfe();
    }

/* ... UNSCALE THE MATRIX, SOLUTION, AND RHS VECTORS. */

L300:
    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___108.ciunit = itcom1_1.nout;
  s_wsfe(&io___108);
  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[5] = itcom3_1.omega;
    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;
} /* sor_ */

/* Subroutine */ int ssorcg_(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  SSORCG\002)";
    static char fmt_40[] = "(\0020\002,\002*** W A R N I N G ************\
\002/\0020\002,\002    IN ITPACK ROUTINE SSORCG\002/\002 \002,\002    RPARM(\
1) =\002,d10.3,\002 (ZETA)\002/\002 \002,\002    A VALUE THIS SMALL MAY HIND\
ER CONVERGENCE \002/\002 \002,\002    SINCE MACHINE PRECISION DRELPR =\002,d\
10.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 SSORCG \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 SSORCG \002/\002\
 \002,\002    ERROR DETECTED IN SUBROUTINE  SBELM \002/\002 \002,\002    WHI\
CH REMOVES ROWS AND COLUMNS OF SYSTEM \002/\002 \002,\002    WHEN DIAGONAL E\
NTRY 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 SSORCG \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 SSORCG \002/\002\
 \002,\002    ERROR DETECTED IN SUBROUTINE  PRBNDX\002/\002 \002,\002    WHI\
CH COMPUTES 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 SSORCG \002/\002\
 \002,\002    ERROR DETECTED IN SUBROUTINE  PERMAT\002/\002 \002,\002    WHI\
CH DOES THE RED-BLACK PERMUTATION\002/\002 \002,\002    IER =