dsrc2c.c

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

    /* 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___1.ciunit = itcom1_1.nout;
  s_wsfe(&io___1);
  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__1);
    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___7.ciunit = itcom1_1.nout;
  s_wsfe(&io___7);
  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 = 11;
    if (itcom1_1.level >= 0) {
  io___12.ciunit = itcom1_1.nout;
  s_wsfe(&io___12);
  do_fio(&c__1, (char *)&n, (ftnlen)sizeof(integer));
  e_wsfe();
    }
    goto L370;
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___14.ciunit = itcom1_1.nout;
  s_wsfe(&io___14);
  do_fio(&c__1, (char *)&ier, (ftnlen)sizeof(integer));
  do_fio(&c__1, (char *)&tol, (ftnlen)sizeof(doublereal));
  e_wsfe();
    }
    goto L370;

/* ... INITIALIZE WKSP BASE ADDRESSES. */

L90:
    ib1 = 1;
    ib2 = ib1 + n;
    ib3 = ib2 + n;
    ib4 = ib3 + n;
    ib5 = ib4 + n;
    iparm[8] = (n << 2) + (itcom1_1.itmax << 1);
    if (itcom1_1.isym != 0) {
  iparm[8] += itcom1_1.itmax << 1;
    }
    if (*nw >= iparm[8]) {
  goto L110;
    }
    ier = 12;
    if (itcom1_1.level >= 0) {
  io___20.ciunit = itcom1_1.nout;
  s_wsfe(&io___20);
  do_fio(&c__1, (char *)&(*nw), (ftnlen)sizeof(integer));
  do_fio(&c__1, (char *)&iparm[8], (ftnlen)sizeof(integer));
  e_wsfe();
    }
    goto L370;

/* ... 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___23.ciunit = itcom1_1.nout;
  s_wsfe(&io___23);
  do_fio(&c__1, (char *)&ier, (ftnlen)sizeof(integer));
  do_fio(&c__1, (char *)&nb, (ftnlen)sizeof(integer));
  e_wsfe();
    }
    goto L370;

/* ... PERMUTE MATRIX AND RHS */

L130:
    if (itcom1_1.level >= 2) {
  io___24.ciunit = itcom1_1.nout;
  s_wsfe(&io___24);
  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___25.ciunit = itcom1_1.nout;
  s_wsfe(&io___25);
  do_fio(&c__1, (char *)&ier, (ftnlen)sizeof(integer));
  e_wsfe();
    }
    goto L370;
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___26.ciunit = itcom1_1.nout;
  s_wsfe(&io___26);
  do_fio(&c__1, (char *)&ier, (ftnlen)sizeof(integer));
  e_wsfe();
    }
    goto L370;
L190:
    if (itcom1_1.level <= 2) {
  goto L220;
    }
    io___27.ciunit = itcom1_1.nout;
    s_wsfe(&io___27);
    e_wsfe();
    if (itcom2_1.adapt) {
  io___28.ciunit = itcom1_1.nout;
  s_wsfe(&io___28);
  e_wsfe();
    }
L220:
    if (iparm[11] != 0) {
  goto L230;
    }
    timi1 = timer_(&dummy);

/* ... COMPUTE INITIAL PSEUDO-RESIDUAL */

L230:
    itpackdcopy_(&n, &rhs[1], &c__1, &wksp[ib2], &c__1);
    pjac_(&n, &ia[1], &ja[1], &a[1], &u[1], &wksp[ib2]);
    vevmw_(&n, &wksp[ib2], &u[1]);

/* ... ITERATION SEQUENCE */

    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 L240;
  }

/* ... CODE FOR THE EVEN ITERATIONS. */

/*     U           = U(IN)             WKSP(IB2) = DEL(IN) */
/*     WKSP(IB1)   = U(IN-1)           WKSP(IB3) = DEL(IN-1) */

  itjcg_(&n, &ia[1], &ja[1], &a[1], &u[1], &wksp[ib1], &wksp[ib2], &
    wksp[ib3], &wksp[ib4], &wksp[ib5]);

  if (itcom2_1.halt) {
      goto L280;
  }
  goto L250;

/* ... CODE FOR THE ODD ITERATIONS. */

/*     U           = U(IN-1)           WKSP(IB2) = DEL(IN-1) */
/*     WKSP(IB1)   = U(IN)             WKSP(IB3) = DEL(IN) */

L240:
  itjcg_(&n, &ia[1], &ja[1], &a[1], &wksp[ib1], &u[1], &wksp[ib3], &
    wksp[ib2], &wksp[ib4], &wksp[ib5]);

  if (itcom2_1.halt) {
      goto L280;
  }
L250:
  ;
    }

/* ... ITMAX HAS BEEN REACHED */

    if (iparm[11] != 0) {
  goto L260;
    }
    timi2 = timer_(&dummy);
    time1 = (doublereal) (timi2 - timi1);
L260:
    ier = 13;
    if (itcom1_1.level >= 1) {
  io___33.ciunit = itcom1_1.nout;
  s_wsfe(&io___33);
  do_fio(&c__1, (char *)&itcom1_1.itmax, (ftnlen)sizeof(integer));
  e_wsfe();
    }
    if (iparm[3] == 0) {
  rparm[1] = itcom3_1.stptst;
    }
    goto L310;

/* ... METHOD HAS CONVERGED */

L280:
    if (iparm[11] != 0) {
  goto L290;
    }
    timi2 = timer_(&dummy);
    time1 = (doublereal) (timi2 - timi1);
L290:
    if (itcom1_1.level >= 1) {
  io___34.ciunit = itcom1_1.nout;
  s_wsfe(&io___34);
  do_fio(&c__1, (char *)&itcom1_1.in, (ftnlen)sizeof(integer));
  e_wsfe();
    }

/* ... PUT SOLUTION INTO U IF NOT ALREADY THERE. */

L310:
    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 L340;
    }
    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 L330;
    }
    if (itcom1_1.level >= 0) {
  io___36.ciunit = itcom1_1.nout;
  s_wsfe(&io___36);
  do_fio(&c__1, (char *)&ierper, (ftnlen)sizeof(integer));
  e_wsfe();
    }
    if (ier == 0) {
  ier = ierper;
    }
    goto L370;
L330:
    pervec_(&n, &rhs[1], &iwksp[ib2]);
    pervec_(&n, &u[1], &iwksp[ib2]);

/* ... OPTIONAL ERROR ANALYSIS */

L340:
    idgts = iparm[12];
    if (idgts < 0) {
  goto L350;
    }
    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 */

L350:
    iparm[8] -= (itcom1_1.itmax - itcom1_1.in) << 1;
    if (iparm[11] != 0) {
  goto L360;
    }
    timj2 = timer_(&dummy);
    time2 = (doublereal) (timj2 - timj1);
L360:
    if (itcom1_1.isym != 0) {
  iparm[8] -= (itcom1_1.itmax - itcom1_1.in) << 1;
    }
    if (iparm[3] != 0) {
  goto L370;
    }
    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;

L370:
    *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;
} /* jcg_ */

/* Subroutine */ int jsi_(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  JSI\002)";
    static char fmt_40[] = "(\0020\002,\002*** W A R N I N G ************\
\002/\0020\002,\002    IN ITPACK ROUTINE JSI\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 JSI \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 JSI \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 JSI \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 ****\