Actual source code: unary.c

      1: #define PETSCMAT_DLL

      3: /* unary.f -- translated by f2c (version of 25 March 1992  12:58:56).

      5:         This code is protected by the GNU copyright. See the file
      6:      ilut.c in this directory for the full copyright. See below for the Author.
      7: */
     8:  #include petsc.h
      9: /* ----------------------------------------------------------------------- */
     10: static PetscErrorCode SPARSEKIT2rperm(PetscInt *nrow,PetscScalar *a,PetscInt *ja,PetscInt *ia,PetscScalar *ao,PetscInt *jao,PetscInt *iao,PetscInt *perm,PetscInt *job)
     11: {
     12:     /* System generated locals */
     13:     PetscInt i__1,i__2;

     15:     /* Local variables */
     16:     PetscInt i,j,k,ii,ko;
     17:     PetscInt values;

     19: /* -----------------------------------------------------------------------
     20:  */
     21: /* this subroutine permutes the rows of a matrix in CSR format. */
     22: /* rperm  computes B = P A  where P is a permutation matrix. */
     23: /* the permutation P is defined through the array perm: for each j, */
     24: /* perm(j) represents the destination row number of row number j. */
     25: /* Youcef Saad -- recoded Jan 28, 1991. */
     26: /* -----------------------------------------------------------------------
     27:  */
     28: /* on entry: */
     29: /* ---------- */
     30: /* n         = dimension of the matrix */
     31: /* a, ja, ia = input matrix in csr format */
     32: /* perm         = integer array of length nrow containing the permutation arrays
     33: */
     34: /*           for the rows: perm(i) is the destination of row i in the */
     35: /*         permuted matrix. */
     36: /*         ---> a(i,j) in the original matrix becomes a(perm(i),j) */
     37: /*         in the output  matrix. */

     39: /* job        = integer indicating the work to be done: */
     40: /*                 job = 1        permute a, ja, ia into ao, jao, iao */
     41: /*                       (including the copying of real values ao and */
     42: /*                       the array iao). */
     43: /*                 job .ne. 1 :  ignore real values. */
     44: /*                     (in which case arrays a and ao are not needed nor
     45: */
     46: /*                      used). */

     48: /* ------------ */
     49: /* on return: */
     50: /* ------------ */
     51: /* ao, jao, iao = input matrix in a, ja, ia format */
     52: /* note : */
     53: /*        if (job.ne.1)  then the arrays a and ao are not used. */
     54: /* ----------------------------------------------------------------------c
     55:  */
     56: /*           Y. Saad, May  2, 1990                                      c
     57: */
     58: /* ----------------------------------------------------------------------c
     59:  */
     60:     /* Parameter adjustments */
     61:     --perm;
     62:     --iao;
     63:     --jao;
     64:     --ao;
     65:     --ia;
     66:     --ja;
     67:     --a;

     69:     /* Function Body */
     70:     values = *job == 1;

     72: /*     determine pointers for output matix. */

     74:     i__1 = *nrow;
     75:     for (j = 1; j <= i__1; ++j) {
     76:         i = perm[j];
     77:         iao[i + 1] = ia[j + 1] - ia[j];
     78: /* L50: */
     79:     }

     81: /* get pointers from lengths */

     83:     iao[1] = 1;
     84:     i__1 = *nrow;
     85:     for (j = 1; j <= i__1; ++j) {
     86:         iao[j + 1] += iao[j];
     87: /* L51: */
     88:     }

     90: /* copying */

     92:     i__1 = *nrow;
     93:     for (ii = 1; ii <= i__1; ++ii) {

     95: /* old row = ii  -- new row = iperm(ii) -- ko = new pointer */

     97:         ko = iao[perm[ii]];
     98:         i__2 = ia[ii + 1] - 1;
     99:         for (k = ia[ii]; k <= i__2; ++k) {
    100:             jao[ko] = ja[k];
    101:             if (values) {
    102:                 ao[ko] = a[k];
    103:             }
    104:             ++ko;
    105: /* L60: */
    106:         }
    107: /* L100: */
    108:     }

    110:     return 0;
    111: /* ---------end-of-rperm -------------------------------------------------
    112:  */
    113: /* -----------------------------------------------------------------------
    114:  */
    115: } /* rperm_ */

    117: /* ----------------------------------------------------------------------- */
    118: static PetscErrorCode SPARSEKIT2cperm(PetscInt *nrow,PetscScalar * a,PetscInt * ja,PetscInt * ia,PetscScalar * ao,PetscInt * jao,PetscInt * iao,PetscInt * perm,PetscInt * job)
    119: {
    120:     /* System generated locals */
    121:     PetscInt i__1;

    123:     /* Local variables */
    124:     PetscInt i,k,nnz;

    126: /* -----------------------------------------------------------------------
    127:  */
    128: /* this subroutine permutes the columns of a matrix a, ja, ia. */
    129: /* the result is written in the output matrix  ao, jao, iao. */
    130: /* cperm computes B = A P, where  P is a permutation matrix */
    131: /* that maps column j into column perm(j), i.e., on return */
    132: /*      a(i,j) becomes a(i,perm(j)) in new matrix */
    133: /* Y. Saad, May 2, 1990 / modified Jan. 28, 1991. */
    134: /* -----------------------------------------------------------------------
    135:  */
    136: /* on entry: */
    137: /* ---------- */
    138: /* nrow         = row dimension of the matrix */

    140: /* a, ja, ia = input matrix in csr format. */

    142: /* perm        = integer array of length ncol (number of columns of A */
    143: /*         containing the permutation array  the columns: */
    144: /*         a(i,j) in the original matrix becomes a(i,perm(j)) */
    145: /*         in the output matrix. */

    147: /* job        = integer indicating the work to be done: */
    148: /*                 job = 1        permute a, ja, ia into ao, jao, iao */
    149: /*                       (including the copying of real values ao and */
    150: /*                       the array iao). */
    151: /*                 job .ne. 1 :  ignore real values ao and ignore iao. */

    153: /* ------------ */
    154: /* on return: */
    155: /* ------------ */
    156: /* ao, jao, iao = input matrix in a, ja, ia format (array ao not needed)
    157: */

    159: /* Notes: */
    160: /* ------- */
    161: /* 1. if job=1 then ao, iao are not used. */
    162: /* 2. This routine is in place: ja, jao can be the same. */
    163: /* 3. If the matrix is initially sorted (by increasing column number) */
    164: /*    then ao,jao,iao  may not be on return. */

    166: /* ----------------------------------------------------------------------c
    167:  */
    168: /* local parameters: */

    170:     /* Parameter adjustments */
    171:     --perm;
    172:     --iao;
    173:     --jao;
    174:     --ao;
    175:     --ia;
    176:     --ja;
    177:     --a;

    179:     /* Function Body */
    180:     nnz = ia[*nrow + 1] - 1;
    181:     i__1 = nnz;
    182:     for (k = 1; k <= i__1; ++k) {
    183:         jao[k] = perm[ja[k]];
    184: /* L100: */
    185:     }

    187: /*     done with ja array. return if no need to touch values. */

    189:     if (*job != 1) {
    190:         return 0;
    191:     }

    193: /* else get new pointers -- and copy values too. */

    195:     i__1 = *nrow + 1;
    196:     for (i = 1; i <= i__1; ++i) {
    197:         iao[i] = ia[i];
    198: /* L1: */
    199:     }

    201:     i__1 = nnz;
    202:     for (k = 1; k <= i__1; ++k) {
    203:         ao[k] = a[k];
    204: /* L2: */
    205:     }

    207:     return 0;
    208: /* ---------end-of-cperm--------------------------------------------------
    209:  */
    210: /* -----------------------------------------------------------------------
    211:  */
    212: } /* cperm_ */

    214: /* ----------------------------------------------------------------------- */
    215: PetscErrorCode SPARSEKIT2dperm(PetscInt *nrow,PetscScalar *a,PetscInt *ja,PetscInt *ia,PetscScalar *ao,PetscInt *jao,PetscInt *iao,PetscInt *perm,PetscInt *qperm,PetscInt *job)
    216: {
    217:     PetscInt locjob;

    219: /* -----------------------------------------------------------------------
    220:  */
    221: /* This routine permutes the rows and columns of a matrix stored in CSR */

    223: /* format. i.e., it computes P A Q, where P, Q are permutation matrices.
    224: */
    225: /* P maps row i into row perm(i) and Q maps column j into column qperm(j):
    226:  */
    227: /*      a(i,j)    becomes   a(perm(i),qperm(j)) in new matrix */
    228: /* In the particular case where Q is the transpose of P (symmetric */
    229: /* permutation of A) then qperm is not needed. */
    230: /* note that qperm should be of length ncol (number of columns) but this
    231: */
    232: /* is not checked. */
    233: /* -----------------------------------------------------------------------
    234:  */
    235: /* Y. Saad, Sep. 21 1989 / recoded Jan. 28 1991. */
    236: /* -----------------------------------------------------------------------
    237:  */
    238: /* on entry: */
    239: /* ---------- */
    240: /* n         = dimension of the matrix */
    241: /* a, ja, */
    242: /*    ia = input matrix in a, ja, ia format */
    243: /* perm         = integer array of length n containing the permutation arrays */
    244: /*           for the rows: perm(i) is the destination of row i in the */
    245: /*         permuted matrix -- also the destination of column i in case */
    246: /*         permutation is symmetric (job .le. 2) */

    248: /* qperm        = same thing for the columns. This should be provided only */
    249: /*         if job=3 or job=4, i.e., only in the case of a nonsymmetric */
    250: /*           permutation of rows and columns. Otherwise qperm is a dummy */

    252: /* job        = integer indicating the work to be done: */
    253: /* * job = 1,2 permutation is symmetric  Ao :== P * A * transp(P) */
    254: /*                 job = 1        permute a, ja, ia into ao, jao, iao */
    255: /*                 job = 2 permute matrix ignoring real values. */
    256: /* * job = 3,4 permutation is non-symmetric  Ao :== P * A * Q */
    257: /*                 job = 3        permute a, ja, ia into ao, jao, iao */
    258: /*                 job = 4 permute matrix ignoring real values. */

    260: /* on return: */
    261: /* ----------- */
    262: /* ao, jao, iao = input matrix in a, ja, ia format */

    264: /* in case job .eq. 2 or job .eq. 4, a and ao are never referred to */
    265: /* and can be dummy arguments. */
    266: /* Notes: */
    267: /* ------- */
    268: /*  1) algorithm is in place */
    269: /*  2) column indices may not be sorted on return even  though they may be
    270:  */
    271: /*     on entry. */
    272: /* ----------------------------------------------------------------------c
    273:  */
    274: /* local variables */

    276: /*     locjob indicates whether or not real values must be copied. */

    278:     /* Parameter adjustments */
    279:     --qperm;
    280:     --perm;
    281:     --iao;
    282:     --jao;
    283:     --ao;
    284:     --ia;
    285:     --ja;
    286:     --a;

    288:     /* Function Body */
    289:     locjob = *job % 2;

    291: /* permute rows first */

    293:     SPARSEKIT2rperm(nrow, &a[1], &ja[1], &ia[1], &ao[1], &jao[1], &iao[1], &perm[1], &locjob);

    295: /* then permute columns */

    297:     locjob = 0;

    299:     if (*job <= 2) {
    300:         SPARSEKIT2cperm(nrow, &ao[1], &jao[1], &iao[1], &ao[1], &jao[1], &iao[1], &perm[1], &locjob);
    301:     } else {
    302:         SPARSEKIT2cperm(nrow, &ao[1], &jao[1], &iao[1], &ao[1], &jao[1], &iao[1], &qperm[1], &locjob);
    303:     }

    305:     return 0;
    306: /* -------end-of-dperm----------------------------------------------------
    307:  */
    308: /* -----------------------------------------------------------------------
    309:  */
    310: } /* dperm_ */

    312: /* ----------------------------------------------------------------------- */
    313: PetscErrorCode SPARSEKIT2msrcsr(PetscInt *n,PetscScalar * a,PetscInt * ja,PetscScalar * ao,PetscInt * jao,PetscInt * iao,PetscScalar * wk,PetscInt * iwk)
    314: {
    315:     /* System generated locals */
    316:     PetscInt i__1, i__2;

    318:     /* Local variables */
    319:     PetscInt iptr;
    320:     PetscInt added;
    321:     PetscInt i, j, k, idiag, ii;

    323: /* -----------------------------------------------------------------------
    324:  */
    325: /*       Modified - Sparse Row  to   Compressed Sparse Row */

    327: /* -----------------------------------------------------------------------
    328:  */
    329: /* converts a compressed matrix using a separated diagonal */
    330: /* (modified sparse row format) in the Compressed Sparse Row */
    331: /* format. */
    332: /* does not check for zero elements in the diagonal. */

    335: /* on entry : */
    336: /* --------- */
    337: /* n          = row dimension of matrix */
    338: /* a, ja      = sparse matrix in msr sparse storage format */
    339: /*              see routine csrmsr for details on data structure */

    341: /* on return : */
    342: /* ----------- */

    344: /* ao,jao,iao = output matrix in csr format. */

    346: /* work arrays: */
    347: /* ------------ */
    348: /* wk       = real work array of length n */
    349: /* iwk      = integer work array of length n+1 */

    351: /* notes: */
    352: /*   The original version of this was NOT in place, but has */
    353: /*   been modified by adding the vector iwk to be in place. */
    354: /*   The original version had ja instead of iwk everywhere in */
    355: /*   loop 500.  Modified  Sun 29 May 1994 by R. Bramley (Indiana). */

    357: /* -----------------------------------------------------------------------
    358:  */
    359:     /* Parameter adjustments */
    360:     --iwk;
    361:     --wk;
    362:     --iao;
    363:     --jao;
    364:     --ao;
    365:     --ja;
    366:     --a;

    368:     /* Function Body */
    369:     i__1 = *n;
    370:     for (i = 1; i <= i__1; ++i) {
    371:         wk[i] = a[i];
    372:         iwk[i] = ja[i];
    373: /* L1: */
    374:     }
    375:     iwk[*n + 1] = ja[*n + 1];
    376:     iao[1] = 1;
    377:     iptr = 1;
    378: /* --------- */
    379:     i__1 = *n;
    380:     for (ii = 1; ii <= i__1; ++ii) {
    381:         added = 0;
    382:         idiag = iptr + (iwk[ii + 1] - iwk[ii]);
    383:         i__2 = iwk[ii + 1] - 1;
    384:         for (k = iwk[ii]; k <= i__2; ++k) {
    385:             j = ja[k];
    386:             if (j < ii) {
    387:                 ao[iptr] = a[k];
    388:                 jao[iptr] = j;
    389:                 ++iptr;
    390:             } else if (added) {
    391:                 ao[iptr] = a[k];
    392:                 jao[iptr] = j;
    393:                 ++iptr;
    394:             } else {
    395: /* add diag element - only reserve a position for it. */
    396:                 idiag = iptr;
    397:                 ++iptr;
    398:                 added = 1;
    399: /*     then other element */
    400:                 ao[iptr] = a[k];
    401:                 jao[iptr] = j;
    402:                 ++iptr;
    403:             }
    404: /* L100: */
    405:         }
    406:         ao[idiag] = wk[ii];
    407:         jao[idiag] = ii;
    408:         if (! added) {
    409:             ++iptr;
    410:         }
    411:         iao[ii + 1] = iptr;
    412: /* L500: */
    413:     }
    414:     return 0;
    415: } /* msrcsr_ */