matrix_skyline.c

有限元程序

              free((char *) d2);
           }
              ZeroUnits(&(spDisp->spColUnits[0]));

         }
           /* [a] : Forward Substitution */

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

           for(ij = 0; ij < spDisp->iNoColumns; ij++)
               for(ii=0; ii < spDisp->iNoRows; ii++)
                   spDisp->uMatrix.daa[ii][ij] /= spA->uMatrix.daa[ii][1];
   
           /* [b] : Back Substitution */

           for(ij = 0; ij < spDisp->iNoColumns; ij++)
               for(ii = spDisp->iNoRows-1; ii >= 0; ii--)
                   for(ik = ii+1; ik <= spDisp->iNoRows-1; ik++)
                       if((ik+1-ii) <= spA->uMatrix.daa[ik][0] )
                           spDisp->uMatrix.daa[ii][ij] -= spDisp->uMatrix.daa[ik][ij] *
                                                       spA->uMatrix.daa[ik][ik+1-ii];

           break;
        default:
             FatalError("In LUBacksubstitutionSkyline() : Undefined Matrix Data Type",
                       (char *) NULL);
    }

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

    return ( spDisp );
}

/*
 *  ===================================================================
 *  MatrixInverseSkyline() : Compute inverse of Skyline Double Matrix
 * 
 *  Input  :  MATRIX *spA -- Pointer to matrix data structure [A].
 *  Output :  MATRIX *spB -- Pointer to matrix inverse [A]^{-1}.
 *  ===================================================================
 */

#ifdef __STDC__
MATRIX *MatrixInverseSkyline ( MATRIX *spA )
#else  /* Start case not STDC */
MATRIX *MatrixInverseSkyline ( spA )
MATRIX *spA;
#endif /* End case not STDC */
{
MATRIX     *spLU;
MATRIX        *F;
MATRIX   *spAinv;
int       ii, ij;
int UNITS_SWITCH;

     UNITS_SWITCH = CheckUnits();

    /* [a] : Check Input */

       if( spA == NULL )
           FatalError("In MatrixInverseSkyline () : [ma] == NULL \n",
                     (char *) NULL);

    /* [b] : Compute LU Decomposition and column vectors of inverse */

       F = MatrixAllocIndirect( spA->cpMatrixName, DOUBLE_ARRAY,
                                spA->iNoRows, spA->iNoColumns);

       for(ii = 0; ii < spA->iNoRows; ii++) {
           for(ij = 0; ij < spA->iNoColumns; ij++) {
               if(ii == ij)
                  F->uMatrix.daa[ii][ij] = 1.0;
               else 
                  F->uMatrix.daa[ii][ij] = 0.0;
           }
       }

       if( UNITS_SWITCH == ON ) {
          for(ii = 1; ii <= spA->iNoRows; ii++)
              ZeroUnits(&(F->spRowUnits[ii-1]));
          for(ii = 1; ii <= spA->iNoColumns; ii++)
              ZeroUnits(&(F->spColUnits[ii-1]));
       }

       spLU   = LUDecompositionSkyline(spA);
       spAinv = LUBacksubstitutionSkyline( spLU, F );

       spAinv = MatrixIndirectToSkyline(spAinv);
       spAinv = MatrixReallocSkyline(spAinv);

       /* Assign units to spAinv */

       if( UNITS_SWITCH == ON ) {
          for(ii = 1; ii <= spA->iNoRows; ii++){
              UnitsPowerRep( &(spAinv->spColUnits[ii-1]), &(spA->spRowUnits[ii-1]), -1.0, YES );
              UnitsPowerRep( &(spAinv->spRowUnits[ii-1]), &(spA->spColUnits[ii-1]), -1.0, YES );
          }
       }

       MatrixFreeSkyline(spLU);
       MatrixFreeIndirect(F);

       return (spAinv);
} 


/*
 *  =================================================================
 *  Subspace() : Solve [A][x] = [Lambda][B][x] by Subspace Iteration.
 *             : Assume [A] and [B] are stored in Skyline Form, and 
 *             : that we will compute "p = iNoVectors" eigenvectors.
 * 
 *             : [A] and [B] are large (nxn) matrices
 *             : Size of [X_{k}] = [X_{k+1}]   = (nxm) matrix
 *             : Size of [Y_{k}] = [Y_{k+1}]   = (nxm) matrix
 *             : Size of [A_{k+1}] = [B_{k+1}] = (mxm) matrix
 *             : Size of [Eigenvalue_{k+1}] = (mxm) diagonal matrix
 *             : Size of [Q_{k+1}]          = (mxm) matrix
 * 
 *  Note : [A] and [B] are stored in skyline form. All other working  
 *         matrices have INDIRECT storage patterns.
 * 
 *  Input :  spEigenvalue  = Allocated Vector for Eigenvalues
 *           spEigenvector = Allocated Vector for Eigenvectors
 *           iNoVectors    = No Eigenvectors to be computed. 
 *  Output : spEigenvalue  = Allocated Vector for Eigenvalues
 *           spEigenvector = Allocated Vector for Eigenvectors
 * 
 *  Written By : M. Austin                             November 1993. 
 *  =================================================================
 */

enum { MAX_SUBSPACE_ITERATIONS = 50 };

#define TOLERANCE  0.01

Subspace( spA, spB, spEigenvalue, spEigenvector , iNoEigen)
MATRIX   *spA,    *spB;
MATRIX   *spEigenvalue;
MATRIX  *spEigenvector;
int           iNoEigen;
{
MATRIX    *spEigenvectorQ;
MATRIX *spAwork, *spBwork;
MATRIX           *spTemp1;
MATRIX            *spA_LU;
MATRIX               *spY;
MATRIX            *spYnew;
MATRIX   *spEigenvalueOld;
MATRIX               *spV;
MATRIX              *spVk;
int     iSize, ii, ij, ik;
int   iMin, iMax,   iLoop;
int  iSubspaceConvergence;
double dEigenvalueOld;
double dConvergence;
double dEigenvalue;
double dMaxValue;

       printf("\n*** Enter Subspace() : Compute %3d Eigenvalues and Eigenvectors\n",
              iNoEigen);

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

       if( spA == NULL || spB == NULL )
           FatalError("In Subpace() : Pointer spA == NULL or spB == NULL",
                     (char *) NULL);

       if( spA->eType != DOUBLE_ARRAY || spB->eType != DOUBLE_ARRAY ) 
           FatalError("In Subspace() : spA->eType != DOUBLE_ARRAY or spB->eType != DOUBLE_ARRAY",
                     (char *) NULL);

       if( spA->eRep != SKYLINE || spB->eRep != SKYLINE ) 
           FatalError("In Subspace() : spA->eRep != SKYLINE or spB->eRep != SKYLINE",
                     (char *) NULL);

    /* [b] : Initialize Starting Eigenvectors [X_{k}] */

       iSize    = spA->iNoRows;
       for(ii = 1; ii <= iNoEigen; ii = ii + 1)
           spEigenvector->uMatrix.daa[ii-1][ii-1] = 1.0;

    /* [c] : Loops for Subspace Iteration */

       spAwork        = MatrixAllocIndirect("[A_{k+1}]", DOUBLE_ARRAY,
                                             iNoEigen, iNoEigen);
       spBwork        = MatrixAllocIndirect("[B_{k+1}]", DOUBLE_ARRAY,
                                             iNoEigen, iNoEigen);

       spEigenvectorQ = MatrixAllocIndirect("Eigenvector [Q]", DOUBLE_ARRAY,
                                             iNoEigen, iNoEigen);

       spVk   = MatrixAllocIndirect("[V_{k}]",   DOUBLE_ARRAY, iSize, iNoEigen);
       spY    = MatrixAllocIndirect("[Y_{k}]",   DOUBLE_ARRAY, iSize, iNoEigen);
       spYnew = MatrixAllocIndirect("[Y_{k+1}]", DOUBLE_ARRAY, iSize, iNoEigen);

       spEigenvalueOld = MatrixCopy( spEigenvalue );

       iLoop = 0; iSubspaceConvergence = FALSE;
       while( iLoop <= (int) MAX_SUBSPACE_ITERATIONS &&
              iSubspaceConvergence == FALSE ) {
           iLoop = iLoop + 1;

           printf("*** In Subspace() : Start Iteration %3d \n", iLoop);

           /* [c.1] : Compute [Y_{k}] = [B]*[X_{k}] */

           for(ii = 1; ii <= iSize; ii = ii + 1) 
           for(ij = 1; ij <= iNoEigen; ij = ij + 1) {
               spY->uMatrix.daa[ ii-1 ][ ij-1 ] = 0.0;
               for(ik = 1; ik <= iSize; ik = ik + 1)  {

                   iMin = (int) MIN( ii, ik );
                   iMax = (int) MAX( ii, ik );

                   if((iMax-iMin+1) <= spB->uMatrix.daa[iMax-1][0]) {
                       spY->uMatrix.daa[ ii-1 ][ ij-1 ] += 
                            spB->uMatrix.daa[ iMax-1 ][ iMax-iMin+1 ] *
                            spEigenvector->uMatrix.daa[ ik-1 ][ ij-1 ] ; 
                   }
               }
           }

           /* [c.2] : Solve [A][V_{k+1}] = [Y_{k}] */

           if(iLoop == 1) {
              spA_LU  = LUDecompositionSkyline(spA);
           }

           for(ii = 1; ii <= iSize; ii = ii + 1) 
           for(ij = 1; ij <= iNoEigen; ij = ij + 1)
               spVk->uMatrix.daa[ ii-1 ][ ij-1 ] = spY->uMatrix.daa[ ii-1 ][ ij-1 ];

           if(iLoop != 1) {
              MatrixFree (spTemp1);
              MatrixFree (spV);
           }
           spV = LUBacksubstitutionSkyline( spA_LU , spVk );

           /* [c.3] : Compute [Y_{new}] = [A] . [V_{k+1}] */

           for(ii = 1; ii <= iSize; ii = ii + 1) 
           for(ij = 1; ij <= iNoEigen; ij = ij + 1) {
               spYnew->uMatrix.daa[ ii-1 ][ ij-1 ] = 0; 
               for(ik = 1; ik <= iSize; ik = ik + 1)  {

                   iMin = (int) MIN( ii, ik );
                   iMax = (int) MAX( ii, ik );

                   if((iMax-iMin+1) <= spA->uMatrix.daa[iMax-1][0]) {
                       spYnew->uMatrix.daa[ ii-1 ][ ij-1 ] += 
                               spA->uMatrix.daa[ iMax-1 ][ iMax-iMin+1 ] *
                               spV->uMatrix.daa[ ik-1 ][ ij-1 ] ; 
                   }
               }
           }

           /* [c.4] : Compute [A_{k+1}] = [V_{k+1}]^T . [Y_{new}] */

           for(ii = 1; ii <= iNoEigen; ii = ii + 1) 
           for(ij = 1; ij <= iNoEigen; ij = ij + 1) {
               spAwork->uMatrix.daa[ ii-1 ][ ij-1 ] = 0.0; 
               for(ik = 1; ik <= iSize; ik = ik + 1) 
                   spAwork->uMatrix.daa[ ii-1 ][ ij-1 ] += 
                            spV->uMatrix.daa[ ik-1 ][ ii-1 ] * 
                            spYnew->uMatrix.daa[ ik-1 ][ ij-1 ];
           }

           /* [c.5] : Compute [Y_{new}] = [B] . [V_{k+1}] */

           for(ii = 1; ii <= iSize; ii = ii + 1) 
           for(ij = 1; ij <= iNoEigen; ij = ij + 1) {
               spYnew->uMatrix.daa[ ii-1 ][ ij-1 ] = 0; 
               for(ik = 1; ik <= iSize; ik = ik + 1)  {

                   iMin = (int) MIN( ii, ik );
                   iMax = (int) MAX( ii, ik );

                   if((iMax-iMin+1) <= spB->uMatrix.daa[iMax-1][0]) {
                       spYnew->uMatrix.daa[ ii-1 ][ ij-1 ] += 
                               spB->uMatrix.daa[ iMax-1 ][ iMax-iMin+1 ] *
                               spV->uMatrix.daa[ ik-1 ][ ij-1 ] ; 
                   }
               }
           }

           /* [c.6] : Compute [B_{k+1}] = [V_{k+1}]^T . [Y_{k+1}] */

           for(ii = 1; ii <= iNoEigen; ii = ii + 1) 
           for(ij = 1; ij <= iNoEigen; ij = ij + 1) {
               spBwork->uMatrix.daa[ ii-1 ][ ij-1 ] = 0.0; 
               for(ik = 1; ik <= iSize; ik = ik + 1) 
                   spBwork->uMatrix.daa[ ii-1 ][ ij-1 ] += 
                            spV->uMatrix.daa[ ik-1 ][ ii-1 ] * 
                            spYnew->uMatrix.daa[ ik-1 ][ ij-1 ];
           }

           /* [c.7] : Solve [A_{k+1}].[Q_{k+1}] = [B_{k+1}].[Q_{k+1}].[Lambda_{k+1}] */

           GeneralizedEigen( spAwork, spBwork, spEigenvalue, spEigenvectorQ );

           /* [c.8] : Check Convergence Criteria */

           iSubspaceConvergence = TRUE;
           for( ii = 1; ii <= spEigenvalue->iNoRows; ii = ii + 1) {
                dEigenvalue    = spEigenvalue->uMatrix.daa[ii-1][0];
                dEigenvalueOld = spEigenvalueOld->uMatrix.daa[ii-1][0];

                dConvergence = ABS(dEigenvalueOld - dEigenvalue)/dEigenvalue;
                if(dConvergence > TOLERANCE) {
                   iSubspaceConvergence = FALSE;
                }
           }

           MatrixFree ( spEigenvalueOld );
           spEigenvalueOld = MatrixCopy( spEigenvalue );

           /* [c.9] : Update [X_{k+1}] = [V_{k+1}].[Q_{k+1}] */

              spTemp1 = MatrixMult( spV, spEigenvectorQ );

           /* [c.10] : Scale [X_{k+1}] so that Max Value is 1 */

              for(ij = 1; ij <= iNoEigen; ij = ij + 1)  {
                  dMaxValue = 0.0;
                  for(ii = 1; ii <= iSize; ii = ii + 1) 
                      if(ABS(spTemp1->uMatrix.daa[ ii-1 ][ ij-1 ]) >= ABS(dMaxValue))
                         dMaxValue = spTemp1->uMatrix.daa[ ii-1 ][ ij-1 ];

                  for(ii = 1; ii <= iSize; ii = ii + 1) 
                      spEigenvector->uMatrix.daa[ ii-1 ][ ij-1 ] =
                                     spTemp1->uMatrix.daa[ ii-1 ][ ij-1 ]/dMaxValue;
              }

       }

       /* [d] : Summarize Algorithm Performance */

       printf("\n*** SUBSPACE ITERATION CONVERGED IN %2d ITERATIONS \n", iLoop);

       /* [e] : Clean Up before Leaving */

       MatrixFree ( spTemp1 );
       MatrixFree ( spV );
       MatrixFree ( spA_LU );
       MatrixFree ( spAwork );
       MatrixFree ( spBwork );
       MatrixFree ( spEigenvectorQ );
       MatrixFree ( spEigenvalueOld );

       MatrixFree ( spYnew );
       MatrixFree ( spY );
       MatrixFree ( spVk );
}