matrix_skyline.c

有限元程序

     UNITS_SWITCH = CheckUnits();

    /* [a] : Check input for [ma] */

       if(spA == NULL)
          FatalError("In MatrixCopySkyline() : [A] = NULL",
                    (char *) NULL);

    /* [b] : Copy Matrix [spCopy] = [spA] */

       ld = iVectorAlloc ( spA->iNoColumns); 
       for(ii = 1; ii <= spA->iNoColumns; ii++)
           ld[ii-1] = spA->uMatrix.daa[ii-1][0];

       spCopy = MatrixAllocSkyline( spA->cpMatrixName, DOUBLE_ARRAY,
                                    spA->iNoRows, spA->iNoColumns, ld);

       switch( UNITS_SWITCH == ON ) {
         case ON:
            for(ii = 1; ii <= spA->iNoRows; ii++) {
                UnitsCopy(&(spCopy->spRowUnits[ii-1]), &(spA->spRowUnits[ii-1]));
                for(ij = 1; ij <= spA->iNoColumns; ij++) {
                    if( ij <= ld[ii-1] )
                        spCopy->uMatrix.daa[ii-1][ij] = spA->uMatrix.daa[ii-1][ij];
                }
            }
            for(ij = 1; ij <= spA->iNoColumns; ij++)
                UnitsCopy(&(spCopy->spColUnits[ij-1]), &(spA->spColUnits[ij-1]));
            break;
        case OFF:
            for(ii = 1; ii <= spA->iNoRows; ii++)
                for(ij = 1; ij <= ld[ii-1]; ij++)
                   spCopy->uMatrix.daa[ii-1][ij] = spA->uMatrix.daa[ii-1][ij];
           break;
        default:
           break;
      }

       free((char *) ld);
       return ( spCopy ); 
}


/* 
 *  ===========================================================
 *  Functions to Cast Matrix Storage Patterns 
 * 
 *     [a] : Indirect matrix to skyline matrix.
 *     [b] : Skyline matrix to indirect matrix.
 * 
 *  Input  :  MATRIX *spA  -- Pointer to matrix data structure.
 *  Output :  MATRIX *spB  -- Pointer to recast matrix.
 *  ===========================================================
 */

#ifdef __STDC__
MATRIX *MatrixIndirectToSkyline( MATRIX *spA )
#else  /* Start case not STDC */
MATRIX *MatrixIndirectToSkyline( spA )
MATRIX *spA;
#endif /* End case not STDC */
{
double **dpB;
int      *ld;
int   ii, ij;
static double REALTOL = 1E-7;

    /* [a] : Check input matrix [A] is defined and symmetric */

       if(spA == NULL)
          FatalError("In MatrixIndirectToSkyline() : spA == NULL",
                    (char *) NULL);

       for(ii = 1; ii <= spA->iNoColumns-1; ii++) {
          for(ij = 0; ij <= ii-1; ij++){
            if(ABS( spA->uMatrix.daa[ij][ii] - spA->uMatrix.daa[ii][ij]) > REALTOL) {
                printf("In MatrixIndirectToSkyline() : %s[%d][%d]-%s[%d][%d] = %le\n",
                       spA->cpMatrixName, ij, ii, spA->cpMatrixName,ii, ij, spA->uMatrix.daa[ij][ii] - spA->uMatrix.daa[ii][ij]);
                printf("In MatrixIndirectToSkyline() : [%s] is not symmetric\n", spA->cpMatrixName);
                FatalError("In MatrixIndirectToSkyline() :", (char *) NULL);
            }
          }
       }

    /* [c] : Compute required Ld[] matrix */

       ld = iVectorAlloc( spA->iNoColumns );
       for(ii = 1; ii <= spA->iNoColumns; ii++) {
       for(ij = 1; (ij < ii)&&(ABS(spA->uMatrix.daa[ij-1][ii-1]) < MINREAL); ij++);
           ld[ii-1] = ii-ij+1;
       }

       switch((int) spA->eType) {
          case DOUBLE_ARRAY:
               dpB = (double **)MyCalloc(spA->iNoRows,sizeof(double*));
               for(ii = 1; ii <= spA->iNoRows; ii++ ) {
                   dpB[ii-1] = (double *) MyCalloc((ld[ii-1]+1), sizeof(double));
                   dpB[ii-1][0] = ld[ii-1];
                   for(ij = 1; ij <= ld[ii-1]; ij++)
                       dpB[ii-1][ij] = spA->uMatrix.daa[ii-ij][ii-1];
               }
               for(ii = 1; ii <= spA->iNoRows; ii++)
                   free ((char *) (spA->uMatrix.daa[ ii-1 ]));
               free ((char *) (spA->uMatrix.daa));
               break;
          default:
               FatalError("In MatrixIndirectToSkyline() : Undefined spA->eType", (char *) NULL);
               break;
      }
      spA->uMatrix.daa = dpB;
      spA->eRep = SKYLINE;
      free((char *) ld);

       return ( spA );
}

/*
 *  ===================================================================
 *  MatrixSkylineToIndirect() : Cast Skyline matrix into Indirect Form.
 *  ===================================================================
 */

#ifdef __STDC__
MATRIX *MatrixSkylineToIndirect( MATRIX *spA )
#else  /* Start case not STDC */
MATRIX *MatrixSkylineToIndirect( spA )
MATRIX *spA;
#endif /* End case not STDC */
{
int  ii, ij, ik;
double    **dpB;

    dpB = MatrixAllocIndirectDouble( spA->iNoRows, spA->iNoColumns );
    for( ii = 1; ii <= spA->iNoColumns; ii++) 
        for( ij = 1; ij <= spA->uMatrix.daa[ii][0]; ij++)
             ik = ABS(ii-ij)+1;
             dpB[ii-1][ik-1] = dpB[ik-1][ii-1] = spA->uMatrix.daa[ii-1][ij];

    for( ii = 1; ii <= spA->iNoColumns; ii++) 
        free((char *) spA->uMatrix.daa[ii-1]);
    free(spA->uMatrix.daa);

    spA->uMatrix.daa = dpB;
    spA->eRep = INDIRECT;

    return (spA);
}

/*
 *  =======================================================================
 *  MatrixScaleSkyline() : Multiply Matrix by Scalar
 * 
 *  Input  :  MATRIX *spA        -- Pointer to matrix data structure [A].
 *         :  double scale       -- double : scale factor c.
 *  Output :  MATRIX *spB        -- Pointer to scaled matrix [B] = c.[A].
 *  =======================================================================
 */

#ifdef __STDC__
MATRIX *MatrixScaleSkyline( MATRIX *spA, double dScale)
#else  /* Start case not STDC */
MATRIX *MatrixScaleSkyline( spA, dScale)
MATRIX   *spA;
double dScale;
#endif /* End case not STDC */
{
MATRIX     *spB;
int *ld, ii, ij;

    /* [a] : Check Input matrix [A] */

       if(spA == (MATRIX *) NULL) 
          FatalError("In MatrixScaleToSkyline() : [A] == NULL",
                    (char *) NULL);

    /* [b] : Scale Skyline Matrix */

       spB = MatrixCopySkyline( spA );

       for( ii = 1; ii <= spA->iNoRows; ii++) 
           for( ij = 1; ij <= spA->uMatrix.daa[ii-1][0]; ij++)
                spB->uMatrix.daa[ii-1][ij] = dScale * spA->uMatrix.daa[ii-1][ij];

       spB = MatrixReallocSkyline ( spB);

       return ( spB);
}

#ifdef __STDC__
double MatrixContentScaleSkyline(MATRIX *spA, int row_no, int col_no )
#else
double MatrixContentScaleSkyline(spA, row_no, col_no )
MATRIX  *spA;
int row_no;        /* row number    */
int col_no;        /* column number */
#endif
{
int iMin, iMax;
int     ii, ij;
double      da;

#ifdef DEBUG
       printf("*** Enter MatrixContentScale() : spA->iNoRows    = %4d\n", spA->iNoRows);
       printf("                        : spA->iNoColumns = %4d\n", spA->iNoColumns);
#endif
      if(CheckUnits()==OFF)
         FatalError("You have to set units ON to use this function","In MatrixScale",(char *)NULL );
 
      iMin = MIN( row_no, col_no );
      iMax = MAX( row_no, col_no );
      ii = iMax - 1;
      ij = iMax - iMin + 1;

     /* [a] Scale for Column of matrix spA */

      if( ij <= spA->uMatrix.daa[ii][0] )
            da = spA->uMatrix.daa[ii][ij];
      else
            da = 0.0;

      if(spA->spColUnits[col_no-1].scale_factor != 0) 
         da = da / spA->spColUnits[col_no-1].scale_factor;
      else {
         printf("==> for column %d, scale_factor of %s = 0 \n",
               col_no, spA->spColUnits[col_no-1].units_name);
         FatalError("Fatal error in MatrixContentScaleSkyline(): ",(char *)NULL);
      }

     /* [b] Scale for Row of matrix spA */

      if(spA->spRowUnits[row_no-1].scale_factor != 0) 
         da = da / spA->spRowUnits[row_no-1].scale_factor;
      else {
         printf("==> for Row %d, scale_factor of %s = 0",
                row_no, spA->spColUnits[row_no-1].units_name);
         FatalError("Fatal error in MatrixScaleSkyline(): ",(char *)NULL);
      }
      return (da);

#ifdef DEBUG
       printf("*** Leave MatrixContentScaleSkyline()\n");
#endif
}
                      
/*
 *  =======================================================================
 *  LUDecompositionSkyline() : [L][D][L]^T Factorization of Skyline Matrix
 * 
 *  Input  :  MATRIX *spA   -- Pointer to skyline matrix [A].
 *  Output :  MATRIX *spLU  -- Pointer to decomposed matrix [L][D][L]^T.
 *  =======================================================================
 */

#ifdef __STDC__
MATRIX *LUDecompositionSkyline( MATRIX *spA )
#else  /* Start case not STDC */
MATRIX *LUDecompositionSkyline( spA )
MATRIX *spA;
#endif /* End case not STDC */
{
MATRIX   *spLU;
int ii, ij, ik;
int    min_col;

    /* [a] : Check Input matrix [A] */

       if(spA == NULL)
          FatalError("In LUDecompositionSkyline() : Pointer to matrix is NULL",
                    (char *) NULL);

    /* [b] : Compute [L][D][L]^T decomposition */

       spLU = MatrixCopySkyline( spA ); 
       switch((int) spLU->eType ) {
           case DOUBLE_ARRAY:
                for(ii = 0; ii < spLU->iNoColumns; ii++) {
                    for(ij = spLU->uMatrix.daa[ii][0]; ij > 1; ij = ij - 1) {

                        min_col = (int) MIN((spLU->uMatrix.daa[ii][0]),
                                      (spLU->uMatrix.daa[ii-ij+1][0]+ij-1));        

                        for(ik=ij+1; ik <= min_col; ik++) 
                            spLU->uMatrix.daa[ii][ij] -=
                                spLU->uMatrix.daa[ii][ik] *
                                spLU->uMatrix.daa[ii-ij+1][ik-ij+1] *
                                spLU->uMatrix.daa[ii-ik+1][1];

                        spLU->uMatrix.daa[ii][ij] /= spLU->uMatrix.daa[ii-ij+1][1];
                    }

                    for(ij = spLU->uMatrix.daa[ii][0]; ij > 1; ij--) 
                        spLU->uMatrix.daa[ii][1] -=
                            spLU->uMatrix.daa[ii][ij] *
                            spLU->uMatrix.daa[ii][ij] *
                            spLU->uMatrix.daa[ii-ij+1][1];

                    if(ABS(spLU->uMatrix.daa[ii][1]) <= MINREAL) {
                       printf(" In LUDecompositionSkyline() : %s[%d][1] = %le \n", spLU->cpMatrixName,ii+1,
                                spLU->uMatrix.daa[ii][1]);
                       FatalError("In LUDecompositionSkyline() : Singular Matrix ",
                                 (char *) NULL);
                    }
                }
                break;
           default:
                FatalError("In LUDecompositionSkyline() : Undefined Matrix Data Type",
                          (char *) NULL);
        }
    
    return ( spLU );
}

/*
 *  ==================================================================================
 *  LUBacksubstitutionSkyline() : Forward/Back Substitution of Factored Skyline Matrix
 * 
 *  Input  :  MATRIX *spA -- Pointer to matrix data structure.
 *  Output :  MATRIX *spX -- Pointer to matrix solution
 *  ==================================================================================
 */

#ifdef __STDC__
MATRIX *LUBacksubstitutionSkyline( MATRIX *spA, MATRIX *spB)
#else  /* Start case not STDC */
MATRIX *LUBacksubstitutionSkyline( spA, spB )
MATRIX *spA;
MATRIX *spB;
#endif /* End case not STDC */
{
int            ii, ij, ik;
DIMENSIONS        *d1,*d2;
MATRIX            *spDisp;


#ifdef MYDEBUG
       printf("*** Enter LUBacksubstitutionSkyline()\n");
#endif

    if(spA == NULL) {
       FatalError("In LUBacksubstitutionSkyline() : spA == NULL",
                 (char *) NULL);
    }

    if(spA->iNoColumns != spB->iNoRows) {
       FatalError("In LUBacksubstitutionSkyline() : Dimensions of",
                  "A[][] and b[] are incompatible                ",
                 (char *) NULL);
    }

    spDisp = MatrixCopyIndirectDouble(spB);

    switch((int) spA->eType ) {
        case DOUBLE_ARRAY:

           /* [0] : Calculate the units of displacements */

         if( CheckUnits() == ON ) {
           for (ii = 1; ii <= spDisp->iNoRows; ii++) {
              d1 = UnitsMult( &(spB->spRowUnits[ii-1]), &(spB->spColUnits[0]) );
              d2 = UnitsMult( &(spA->spRowUnits[ii-1]), &(spA->spColUnits[ii-1]) );
              UnitsDivRep(&(spDisp->spRowUnits[ii-1]), d1, d2, YES); 

              free((char *) d1->units_name);
              free((char *) d1);
              free((char *) d2->units_name);