matrix_indirect.c

有限元程序

    spScale = VectorAlloc("Scale Factor",  DOUBLE_ARRAY, spA->iNoRows );

    for (ii = 1; ii <= spA->iNoRows; ii = ii + 1) {
         spScale->uVector.da[ii-1] = ABS(spA->uMatrix.daa[ii-1][0]);
         for (ij = 2; ij <= spA->iNoColumns; ij = ij + 1) {
              dTemp = ABS(spA->uMatrix.daa[ii-1][ij-1]);
              spScale->uVector.da[ii-1] = (double) MAX( spScale->uVector.da[ii-1], dTemp );
         }

         if( spScale->uVector.da[ii-1] == 0)
             FatalError("Execution halted in SetupScaleFactors()",
                        "Matrix A is Singular !!",
                        (char *) NULL);
    }

    return ( spScale );
}

/* 
 *  ==================================================================
 *  SetupPivotVector() : Initialize (nx1) vector for row permutations.
 *  
 *  Input :   MATRIX spA       -- Pointer to (nxn) matrix [A].
 *  Output :  VECTOR spPivot   -- Pointer to (nx1) vector spPivot.
 *  ==================================================================
 */ 

#ifdef __STDC__
VECTOR *SetupPivotVector( MATRIX *spA )
#else  /* Start case not STDC */
VECTOR *SetupPivotVector( spA )
MATRIX *spA;
#endif /* End case not STDC */
{
VECTOR *spPivot;
int          ii;

    spPivot = VectorAlloc(  "Row Pivots", INTEGER_ARRAY, spA->iNoRows );

    for (ii = 1; ii <= spA->iNoRows; ii = ii + 1)
         spPivot->uVector.ia[ii-1] = ii;

    return ( spPivot );
}


/* 
 *  ==================================================================
 *  Pivot() : Find largest element below diagonal and pivot rows
 *  
 *  Input  :   MATRIX spA       -- Pointer to (nxn) matrix [A].
 *         :   VECTOR spScale   -- Pointer to (nx1) vector spScale.
 *         :   VECTOR spPivot   -- Pointer to (nx1) vector spPivot.
 *         :   int iColumnNo    -- Column No for Pivoting
 *  Output :   MATRIX spA       -- [A] and [spPivot] are rearranged by
 *                                 the pivoting procedure.
 *  ==================================================================
 */ 

#ifdef __STDC__
Pivot( MATRIX *spA, VECTOR *spScale, VECTOR *spPivot, int iColumnNo )
#else  /* Start case not STDC */
Pivot( spA, spScale, spPivot, iColumnNo )
VECTOR *spScale, *spPivot;
MATRIX     *spA;
int   iColumnNo;
#endif /* End case not STDC */
{
double dLargest,     dTemp;
int  ii , iTemp,    iOrder;
int      iPivot;

    /* [a] : Find Row No of Largest Scaled Equation  */

       iPivot  = iColumnNo;
       iOrder  = spPivot->uVector.ia[ iColumnNo - 1 ];

       dLargest = spA->uMatrix.daa[ iOrder-1 ][ iColumnNo - 1 ] /
                  spScale->uVector.da[ iOrder-1 ];
       dLargest = ABS( dLargest );

       for (ii = iColumnNo + 1; ii <= spA->iNoColumns; ii = ii + 1) {
            iOrder = spPivot->uVector.ia[ ii - 1 ];

            dTemp  = spA->uMatrix.daa[ iOrder-1 ][ iColumnNo - 1 ] /
                     spScale->uVector.da[ iOrder-1 ];
            dTemp  = ABS( dTemp );

            if( dTemp > dLargest ) {
                dLargest = dTemp;
                iPivot   = ii;
            }
       }

       if( dLargest == 0 )
           FatalError("Execution halted in Pivot()",
                      "Matrix A is Singular !!",
                      (char *) NULL);

    /* [b] : Swap Row Nos in Order Vector */

       iTemp  = spPivot->uVector.ia[ iPivot - 1 ];
       spPivot->uVector.ia[ iPivot - 1 ]    = spPivot->uVector.ia[ iColumnNo - 1 ];
       spPivot->uVector.ia[ iColumnNo - 1 ] = iTemp;
}


/*
 *  =================================================================
 *  GeneralizedEigen() : Solve Eigen Problem [A][x] = [B][x][Lambda].
 *                       where [A] and [B] are small fully populated  
 *                       symmetric matrices.  
 *  
 *  We convert [A][x] = [B][x][Lambda] into [A*][y] = [y][Lambda] in
 *  a four-step procedure:  
 *  
 *        (a) : Compute [B] = [L].[L]^T  
 *        (b) : Solve   [L][C] = [A].
 *        (c) : Compute [C]^T.
 *        (d) : Solve   [L][A*]^T = [C]^T; matrix [A*] is symmetric.  
 *  
 *  Use Householder method to solve [A*][y] = [y][Lamda], followed by 
 *  backward substitutution for [L]^T.[x] = [y].
 *  
 *  Written By : M. Austin                              November 1993 
 *  =================================================================
 */

GeneralizedEigen( spA, spB, spEigenvalue, spEigenvector )
MATRIX   *spA,    *spB;
MATRIX   *spEigenvalue;
MATRIX  *spEigenvector;
{
MATRIX *spEigenvectorY;
MATRIX       *spCwork2;
MATRIX        *spAwork;
MATRIX        *spCwork;
MATRIX          *spLLT;
int        iSize;
int   ii, ij, ik;
double       sum;

#ifdef MYDEBUG
       printf("\n*** Enter GeneralizedEigen()\n");
       MatrixPrint ( spA );
       MatrixPrint ( spB );
#endif

    /* [a] : Check Input and Allocate Working Arrays */

       iSize = spA->iNoRows;
       spCwork = MatrixAllocIndirect("Working [C]", DOUBLE_ARRAY, iSize, iSize );
       spAwork = MatrixAllocIndirect("Working [A]", DOUBLE_ARRAY, iSize, iSize );

    /* [b] : Compute Cholesky Decomposition [B] = [L].[L]^T */

       spLLT = MatrixCopy( spB );
       spLLT = CholeskyDecompositionIndirect( spLLT );

    /* [c] : Solve [L][C] = [A] via Forward Substitution */

       for (ij = 1; ij <= spA->iNoRows; ij = ij + 1) {
            for (ii = 1; ii <= spA->iNoColumns; ii = ii + 1) {

            sum  = 0.0;
            for (ik = 1; ik < ii; ik = ik + 1)
                 sum += spLLT->uMatrix.daa[ ii-1 ][ ik-1 ] *
                               spCwork->uMatrix.daa[ ik-1 ][ij-1];

            spCwork->uMatrix.daa[ii-1][ij-1] = (spA->uMatrix.daa[ ii-1 ][ ij-1 ] - sum) /
                                                spLLT->uMatrix.daa[ ii-1 ][ ii-1 ];

            }
       }

    /* [d] : Replace [C] by [C]^T. */

       spCwork2 = MatrixTranspose ( spCwork );

    /* [e] : Solve [L][A*] = [C]^T via Forward Substitution */

       for (ij = 1; ij <= spA->iNoRows; ij = ij + 1) {
            for (ii = 1; ii <= spA->iNoColumns; ii = ii + 1) {

            sum  = 0.0;
            for (ik = 1; ik < ii; ik = ik + 1)
                 sum += spLLT->uMatrix.daa[ ii-1 ][ ik-1 ] *
                               spAwork->uMatrix.daa[ ik-1 ][ij-1];

            spAwork->uMatrix.daa[ii-1][ij-1] =
                    (spCwork2->uMatrix.daa[ ii-1 ][ ij-1 ] - sum) /
                     spLLT->uMatrix.daa[ ii-1 ][ ii-1 ];
            }
       }

    /* [f] : Use Householder Method to Solve Standard Eigenvalue Problem */

       iSize = spAwork->iNoRows;
       spEigenvectorY = MatrixAllocIndirect("Eigenvector [Y]", DOUBLE_ARRAY, iSize, iSize);

    /* [e] : Use Householder Transformation to convert [spStiff] to Triangular form */

       Householder( spAwork, spEigenvalue, spEigenvectorY );

    /* [g] : Compute [L]^T.[X] = [Y] via Back Substitution */

       for (ij = spAwork->iNoRows; ij >= 1 ; ij = ij - 1) {
            for (ii = 1; ii <= spAwork->iNoColumns; ii = ii + 1) {

            sum  = 0.0;
            for (ik = ij + 1; ik <= spAwork->iNoRows; ik = ik + 1)
                 sum += spLLT->uMatrix.daa[ ij-1 ][ ik-1 ] *
                               spEigenvector->uMatrix.daa[ ik-1 ][ii-1];

            spEigenvector->uMatrix.daa[ij-1][ii-1] =
                          (spEigenvectorY->uMatrix.daa[ ij-1 ][ ii-1 ] - sum) /
                           spLLT->uMatrix.daa[ ij-1 ][ ij-1 ];

            }
       }

    /* [h] : Cleanup before leaving */

       MatrixFree( spAwork );
       MatrixFree( spCwork );
       MatrixFree( spCwork2 );
       MatrixFree( spEigenvectorY );
       MatrixFree( spLLT );

#ifdef MYDEBUG
       printf("\n*** Leave GeneralizedEigen()\n");
#endif

}


/*
 *  ========================================================================
 *  Householder() : Compute Eigenvalues/Eigenvectors of Symmetric Matrix Via
 *  Householder's Transformation and QL-Algorithm.
 *  
 *  Input  :  MATRIX *spK           -- Pointer to large matrix [K].
 *         :  MATRIX *spEigenvalue  -- Pointer to (nx1) eigenvalue matrix.
 *         :  MATRIX *spEigenvector -- Pointer to (nxn) eigenvector matrix.
 *  Output :  MATRIX *spEigenvalue  -- Pointer to (nx1) eigenvalue matrix.
 *         :  MATRIX *spEigenvector -- Pointer to (nxn) eigenvector matrix.
 *  
 *  Written By: Mark Austin                                    November 1993
 *  ========================================================================
 */

Householder( spK, spEigenvalue, spEigenvector )
MATRIX            *spK;
MATRIX   *spEigenvalue;
MATRIX  *spEigenvector;
{
double  *t1, *t2;
int         stat;
int       ii, ij;
int iNoEquations;

 /* [a] : Setup Working Matrices for Eigenvalues and Eigenvectors */

    iNoEquations = spK->iNoRows;
    t1 = (double *) MyCalloc( iNoEquations, sizeof(double));
    t2 = (double *) MyCalloc( iNoEquations, sizeof(double));

    for( ii = 1; ii <= iNoEquations; ii = ii + 1)
    for( ij = 1; ij <= iNoEquations; ij = ij + 1)
         spEigenvector->uMatrix.daa[ ii-1 ][ ij-1 ] = spK->uMatrix.daa[ ii-1 ][ ij-1 ];

 /* [b] : Convert [K] to Tridiagonal Form with Householder Transformation */

    HouseHolderTransformation( spEigenvector, t1, t2, iNoEquations);

 /* [c] : Use ql Algorithm to compute Eigenvalues and Eigenvectors */

    stat = QlAlgorithm( spEigenvector, t1, t2, iNoEquations);

 /* [d] : Transfer Working Matrix to Eigenvalues */

    for( ii = 1; ii <= iNoEquations; ii = ii + 1)
         spEigenvalue->uMatrix.daa[ ii-1 ][0] = t1[ ii-1 ];

    free((char *) t1);
    free((char *) t2);

    return stat;
}


/*
 *  ==============================================================
 *  HouseHolderTrans() : Use HouseHolder Transformation to convert
 *                       Symmetric Matrix to tridiagonal form.
 *  ==============================================================
 */

static double TOLERANCE = 1.3e-16;

HouseHolderTransformation( spK , d, e, iNoEquations)
MATRIX      *spK;
double    *d, *e;
int iNoEquations;
{
int    ii, ij, ik, im;
double dTempf, dTempg;
double dTemph, dTempk;

    ii = iNoEquations;
    while (--ii) {

        im = ii - 2;
        dTempf = spK->uMatrix.daa[ii][ii-1];
        dTempg = 0.0;
        for (ij = 0; ij <= im; ij++)
             dTempg += spK->uMatrix.daa[ii][ij] * spK->uMatrix.daa[ii][ij];

        dTemph = dTempg + dTempf * dTempf;
        if (dTempg < TOLERANCE) {
                 e[ii]  = dTempf;
                 dTemph = 0.0;
        } else {
             im++;
             if (dTempf >= 0)
                 dTempg = -sqrt( dTemph );
             else
                 dTempg =  sqrt( dTemph );

             e[ii] = dTempg;
             dTemph -= dTempf * dTempg;
             spK->uMatrix.daa[ii][ii-1] = dTempf - dTempg;

             dTempf = 0;
             for (ij = 0; ij <= im; ij++) {
                  spK->uMatrix.daa[ij][ii] = spK->uMatrix.daa[ii][ij] / dTemph;
                  dTempg = 0;
                  for (ik = 0; ik <= ij; ik++)
                       dTempg += spK->uMatrix.daa[ij][ik] * spK->uMatrix.daa[ii][ik];
                  for (ik = ij+1; ik <= im; ik++)
                       dTempg += spK->uMatrix.daa[ik][ij] * spK->uMatrix.daa[ii][ik];
                  e[ij] = dTempg / dTemph;
                  dTempf += dTempg * spK->uMatrix.daa[ij][ii];
             }

             dTempk = dTempf / (dTemph + dTemph);
             for (ij = 0; ij <= im; ij++) {
                  dTempf = spK->uMatrix.daa[ii][ij];
                  dTempg = e[ij] - dTempk * dTempf;
                  e[ij]  = dTempg;
                  for (ik=0; ik <= ij; ik++)
                       spK->uMatrix.daa[ij][ik] = spK->uMatrix.daa[ij][ik] - dTempf * e[ik] -
                                                  dTempg * spK->uMatrix.daa[ii][ik];
             }
        }

        d[ii] = dTemph;

    }

    d[0] = 0.0;
    e[0] = 0.0;
    for (ii=0; ii < iNoEquations; ii++) {
         im = ii - 1;
         if (d[ii] != 0.0) 
         for (ij=0; ij <= im; ij++) {
              dTempg = 0.0;
              for (ik=0; ik <= im; ik++)
                   dTempg += spK->uMatrix.daa[ii][ik] * spK->uMatrix.daa[ik][ij];
              for (ik=0; ik <= im; ik++)
                   spK->uMatrix.daa[ik][ij] -= dTempg * spK->uMatrix.daa[ik][ii];
         }

         d[ii] = spK->uMatrix.daa[ii][ii];
         spK->uMatrix.daa[ii][ii] = 1.0;
         for (ij=0; ij <= im; ij++)
              spK->uMatrix.daa[ii][ij] = spK->uMatrix.daa[ij][ii] = 0.0;
    }
}


/* 
 *  ================================================
 *  QL() : QL Algorithm : Adapted from EISPACK......
 *  ================================================
 */

QlAlgorithm( spK, d, e, iNoEquations )
MATRIX      *spK;
doub