dsrc2c.c

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

/**
 * \file itpack.cxx
 * \brief Routines from ITPACK numerical library.
 *
 * In this file routines from ITPACK are defined in the itk::fem::itpack
 * namespace. This is done by calling the original ITPACK code
 * (file dsrc2c.c) which was converted from Fortran using f2c.
 * See "ftp://netlib.bell-labs.com/netlib/f2c" for more info.
 * 
 * Besides that we also define some functions that the ITPACK uses and should
 * normally be defined in f2c.lib. But since we don't have that library
 * we redefine these function here, so that ITPACK code is happy.
 *
 * ITPACK was converted with the following command line options:
 *    f2c -C++P dsrc2c.f
 *
 * This produced two files:
 *   - dsrc2c.c was converted to itpack_dsrc2c.c by removing the
 *     #include "f2c" and all "#ifdef __cplusplus extern C  #endif" parts.
 *   - dsrc2c.P was used to obtain function prototypes which are stored
 *     in itpack.h file.
 *
 * To use ITPACK in your ITK code, simply include the header "itpack.h",
 * and link to the corresponding library.
 *
 * \note All ITPACK functions reside in namespace itpack.
 */




/* Turn off warnings in f2c generated code */
#if defined(_MSC_VER)
#if defined(__ICL)
#pragma warning(disable: 239 264 1011 )
#else
#pragma warning(disable: 4101 4244 4554 4756 4723)
#endif
#endif





#include "f2c.h"




integer i_sign(integer *a, integer *b)
{
integer x;
x = (*a >= 0 ? *a : - *a);
return( *b >= 0 ? x : -x);
}


/*
 * The following couple of functions have something to
 * do with I/O. Since we don't want ITPACK to output
 * anything, the functions do nothing.
 */
integer do_fio(integer * a, char * b, ftnlen c )
{
  (void)a; (void)b; (void)c;
  return 0;
}

integer e_wsfe(void)
{
  return 0;
}

integer s_wsfe(cilist * list)
{
  (void)list;
  return 0;
}

doublereal etime_(float *tarray)
{
  tarray[0]=0.0;
  tarray[1]=0.0;
  return 0.0;
}





/*****   END OF ITK ADDITIONS ****/



/* dsrc2c.f -- translated by f2c (version 20020621).
   You must link the resulting object file with the libraries:
  -lf2c -lm   (in that order)
*/

#ifdef __cplusplus
extern "C" {
#endif
#include "f2c.h"

/* Common Block Declarations */

static struct {
    integer in, is, isym, itmax, level, nout;
} itcom1_;

#define itcom1_1 itcom1_

static struct {
    logical adapt, betadt, caseii, halt, partad;
} itcom2_;

#define itcom2_1 itcom2_

static struct {
    doublereal bdelnm, betab, cme, delnnm, delsnm, ff, gamma, omega, qa, qt, 
      rho, rrr, sige, sme, specr, spr, drelpr, stptst, udnm, zeta;
} itcom3_;

#define itcom3_1 itcom3_

/* Table of constant values */

static integer c__1 = 1;
static integer c__0 = 0;
static doublereal c_b21 = 0.;
static integer c__2 = 2;
static integer c__3 = 3;
static integer c__4 = 4;
static doublereal c_b286 = 1.;
static integer c__5 = 5;
static integer c__6 = 6;
static integer c__7 = 7;

/* Subroutine */ int jcg_(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  JCG\002)";
    static char fmt_40[] = "(\0020\002,\002*** W A R N I N G ************\
\002/\0020\002,\002    IN ITPACK ROUTINE JCG\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 JCG \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 JCG \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 JCG \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 JCG \002/\002 \002,\
\002    ERROR DETECTED IN SUBROUTINE  PRBNDX\002/\002 \002,\002    WHICH COM\
PUTES 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 JCG \002/\002 \002,\
\002    ERROR DETECTED IN SUBROUTINE  PERMAT\002/\002 \002,\002    WHICH DOE\
S THE RED-BLACK PERMUTATION\002/\002 \002,\002    IER = \002,i5)";
    static char fmt_180[] = "(\0020\002,\002*** F A T A L     E R R O R ****\
********\002/\0020\002,\002    CALLED FROM ITPACK ROUTINE JCG \002/\002 \002,\
\002    ERROR DETECTED IN SUBROUTINE  SCAL  \002/\002 \002,\002    WHICH SCA\
LES THE SYSTEM   \002/\002 \002,\002    IER = \002,i5)";
    static char fmt_200[] = "(///1x,\002IN THE FOLLOWING, RHO AND GAMMA AR\
E\002,\002 ACCELERATION PARAMETERS\002)";
    static char fmt_210[] = "(1x,\002CME IS THE ESTIMATE OF THE LARGEST EIGE\
NVALUE OF\002,\002 THE JACOBI MATRIX\002)";
    static char fmt_270[] = "(\0020\002,\002*** W A R N I N G ***********\
*\002/\0020\002,\002    IN ITPACK ROUTINE JCG\002/\002 \002,\002    FAILURE \
TO CONVERGE IN\002,i5,\002 ITERATIONS\002)";
    static char fmt_300[] = "(/1x,\002JCG  HAS CONVERGED IN \002,i5,\002 ITE\
RATIONS\002)";
    static char fmt_320[] = "(\0020\002,\002*** F A T A L     E R R O R ****\
********\002/\0020\002,\002    CALLED FROM ITPACK ROUTINE JCG \002/\002 \002,\
\002    ERROR DETECTED IN SUBROUTINE  PERMAT\002/\002 \002,\002    WHICH UND\
OES THE RED-BLACK PERMUTATION   \002/\002 \002,\002    IER = \002,i5)";

    /* System generated locals */
    integer i__1;

    /* Builtin functions */
    integer s_wsfe(cilist *), e_wsfe(void), do_fio(integer *, char *, ftnlen);

    /* Local variables */
    static integer n, n3, nb, ib1, ib2, ib3, ib4, ib5, ier;
    static doublereal tol;
    extern /* Subroutine */ int pjac_(integer *, integer *, integer *, 
      doublereal *, doublereal *, doublereal *), scal_(integer *, 
      integer *, integer *, doublereal *, doublereal *, doublereal *, 
      doublereal *, integer *, integer *, integer *);
    static doublereal temp;
    static integer loop;
    static doublereal time1, time2;
    static real timi1, timj1, timi2, timj2;
    extern /* Subroutine */ int itjcg_(integer *, integer *, integer *, 
      doublereal *, doublereal *, doublereal *, doublereal *, 
      doublereal *, doublereal *, doublereal *), sbelm_(integer *, 
      integer *, integer *, doublereal *, doublereal *, integer *, 
      doublereal *, doublereal *, integer *, integer *, integer *, 
      integer *);
    static integer idgts;
    extern /* Subroutine */ int vfill_(integer *, doublereal *, doublereal *),
       itpackdcopy_(integer *, doublereal *, integer *, doublereal *, integer 
      *);
    extern doublereal timer_(real *);
    static real dummy;
    extern /* Subroutine */ int vevmw_(integer *, doublereal *, doublereal *);
    static doublereal digit1, digit2;
    static integer itmax1;
    extern /* Subroutine */ int echall_(integer *, integer *, integer *, 
      doublereal *, doublereal *, integer *, doublereal *, integer *), 
      pervec_(integer *, doublereal *, integer *), ivfill_(integer *, 
      integer *, integer *);
    static integer ierper;
    extern /* Subroutine */ int echout_(integer *, doublereal *, integer *), 
      permat_(integer *, integer *, integer *, doublereal *, integer *, 
      integer *, integer *, integer *, integer *, integer *), unscal_(
      integer *, integer *, integer *, doublereal *, doublereal *, 
      doublereal *, doublereal *), prbndx_(integer *, integer *, 
      integer *, integer *, integer *, integer *, integer *, integer *, 
      integer *), perror_(integer *, integer *, integer *, doublereal *,
       doublereal *, doublereal *, doublereal *, doublereal *, 
      doublereal *, integer *);

    /* Fortran I/O blocks */
    static cilist io___1 = { 0, 0, 0, fmt_10, 0 };
    static cilist io___7 = { 0, 0, 0, fmt_40, 0 };
    static cilist io___12 = { 0, 0, 0, fmt_60, 0 };
    static cilist io___14 = { 0, 0, 0, fmt_80, 0 };
    static cilist io___20 = { 0, 0, 0, fmt_100, 0 };
    static cilist io___23 = { 0, 0, 0, fmt_120, 0 };
    static cilist io___24 = { 0, 0, 0, fmt_140, 0 };
    static cilist io___25 = { 0, 0, 0, fmt_150, 0 };
    static cilist io___26 = { 0, 0, 0, fmt_180, 0 };
    static cilist io___27 = { 0, 0, 0, fmt_200, 0 };
    static cilist io___28 = { 0, 0, 0, fmt_210, 0 };
    static cilist io___33 = { 0, 0, 0, fmt_270, 0 };
    static cilist io___34 = { 0, 0, 0, fmt_300, 0 };
    static cilist io___36 = { 0, 0, 0, fmt_320, 0 };



/*     ITPACK 2C MAIN SUBROUTINE  JCG  (JACOBI CONJUGATE GRADIENT) */
/*     EACH OF THE MAIN SUBROUTINES: */
/*           JCG, JSI, SOR, SSORCG, SSORSI, RSCG, RSSI */
/*     CAN BE USED INDEPENDENTLY OF THE OTHERS */

/* ... FUNCTION: */

/*          THIS SUBROUTINE, JCG, DRIVES THE JACOBI CONJUGATE */
/*          GRADIENT ALGORITHM. */

/* ... PARAMETER LIST: */

/*          N      INPUT INTEGER.  DIMENSION OF THE MATRIX. (= NN) */
/*          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.  JACOBI CONJUGATE */
/*                 GRADIENT NEEDS THIS TO BE IN LENGTH AT LEAST */
/*                 4*N + 2*ITMAX,  IF ISYM = 0  (SYMMETRIC STORAGE) */
/*                 4*N + 4*ITMAX,  IF ISYM = 1  (NONSYMMETRIC STORAGE) */
/*                 HERE ITMAX = IPARM(1) AND ISYM = IPARM(5) */
/*                 (ITMAX IS THE MAXIMUM ALLOWABLE NUMBER OF ITERATIONS) */
/*          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) */

/* ... JCG SUBPROGRAM REFERENCES: */

/*          FROM ITPACK    BISRCH, CHGCON, DETERM, DFAULT, ECHALL, */
/*                         ECHOUT, EIGVNS, EIGVSS, EQRT1S, ITERM, TIMER, */
/*                         ITJCG, IVFILL, PARCON, PERMAT, */
/*                         PERROR, PERVEC, PJAC, PMULT, PRBNDX, */
/*                         PSTOP, QSORT, DAXPY, SBELM, SCAL, DCOPY, */
/*                         DDOT, SUM3, UNSCAL, VEVMW, VFILL, VOUT, */
/*                         WEVMW, ZBRENT */
/*          SYSTEM         DABS, DLOG10, DBLE(AMAX0), DMAX1, 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 */
/*                       [phone] */

/*     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 */


/* **** BEGIN: ITPACK COMMON */




/* **** END  : ITPACK COMMON */

/* ... 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 */

/* ... INITIALIZE COMMON BLOCKS */