/* * ====================================================================== * Functions for Operations on Skyline Matrices. * ====================================================================== * * Copyright (C) 1993-96 by Mark Austin and David Mazzoni. * * This software is provided "as is" without express or implied warranty. * Permission is granted to use this software on any computer system, * and to redistribute it freely, subject to the following restrictions: * * 1. The authors are not responsible for the consequences of use of * this software, even if they arise from defects in the software. * 2. The origin of this software must not be misrepresented, either * by explicit claim or by omission. * 3. Altered versions must be plainly marked as such, and must not * be misrepresented as being the original software. * 4. This notice is to remain intact. * * Written by: Lanheng Jin and M.A. Austin July 1992 - October 1993 * ====================================================================== */ #include #include #ifdef __STDC__ #include #include #include #endif #include "defs.h" #include "matrix.h" #include "vector.h" #include "miscellaneous.h" /* #define MYDEBUG */ /* * ======================================================================= * matAllocSkyline() : Allocate (and free) memory for Skyline matrix * data structure. * * Input : char *cpName -- Pointer to name of matrix. * : DATA_TYPE Type -- Data type to be stored in Matrix. * : int iNoRows -- No of Rows in Matrix. * : int iNoColumns -- No of Columns in Matrix. * : int *ld -- Length of each column in skyline profile * Output : MATRIX *spA -- Pointer to matrix data structure. * ======================================================================= */ #ifdef __STDC__ MATRIX *matAllocSkyline( char *cpName, DATA_TYPE Type, int iNoRows, int iNoColumns, int *ld) #else /* start case not STDC */ MATRIX *matAllocSkyline( cpName, Type, iNoRows, iNoColumns, ld) char *cpName; DATA_TYPE Type; int iNoRows; int iNoColumns; int *ld; #endif /* end case not STDC */ { MATRIX *spA; int ii; /* [a] : Check input parameters/arrays */ if(ld == NULL) fatalError("In matAllocSkyline() : ld[] = NULL", (char *) NULL); /* [b] : Allocate matrix parent data structure, and name */ spA = (MATRIX *) safeMalloc( sizeof(MATRIX), __FILE__, __LINE__ ); if(cpName != NULL) spA->cpName = saveString(cpName, __FILE__, __LINE__ ); else spA->cpName = (char *) NULL; /* [c] : Set parameters and allocate memory for matrix array */ spA->Rep = Skyline; spA->Type = Type; spA->iNoRows = iNoRows; spA->iNoColumns = iNoColumns; switch((int) spA->Type) { case DoubleArray: spA->uMatrix.dpp = (double **) safeCalloc( iNoColumns, sizeof(double*), __FILE__, __LINE__ ); for(ii = 0; ii < iNoColumns; ii = ii + 1) { spA->uMatrix.dpp[ ii ] = (double *) safeCalloc((ld[ii]+1), sizeof(double), __FILE__, __LINE__ ); spA->uMatrix.dpp[ ii ][0] = ld[ ii ]; } break; default: fatalError("In matAllocSkyline() : Undefined spA->Type", (char *) NULL); break; } return ( spA ); } /* * ======================================================================= * matFreeSkyline() : Free memory in Skyline Matrix * * Input : MATRIX *Matrix -- Pointer to matrix data structure. * Output : void * ======================================================================= */ #ifdef __STDC__ void matFreeSkyline( MATRIX *spA ) #else /* Start case not STDC */ void matFreeSkyline( spA ) MATRIX *spA; #endif /* End case not STDC */ { int ii; /* [a] : Free body of skyline matrix[][] */ switch( spA->Type ) { case DoubleArray: for(ii = 1; ii <= spA->iNoColumns; ii++) free ((char *) (spA->uMatrix.dpp[ ii ])); free ((char *) (spA->uMatrix.dpp)); break; default: fatalError("In matFreeSkyline() : Undefined spA->Type", (char *) NULL); break; } /* [b] : Free name and parent data structure for matrix */ free ((char *) (spA->cpName)); free ((char *) spA); } /* * ======================================================================= * Print a Matrix [iNoRows][iNoColumns] of data type DOUBLE * Dump Profile of Skyline Matrix. * * Where -- spA->iNoRows = Number of rows in matrix. * spA->iNoColumns = Number of columns in matrix. * COLUMNS_ACROSS_PAGE = Number of columns printed across page. * * ib = Current No of Matrix Block. * iFirstColumn = Number of first column in block * iLastColumn = Number of last column in block * * Input : MATRIX *spA -- Pointer to matrix data structure. * Output : void * ======================================================================= */ enum { COLUMNS_ACROSS_PAGE = 6 }; /* Item [a] */ #ifdef __STDC__ void matPrintSkylineDouble( MATRIX *spA ) #else /* Start case not STDC */ void matPrintSkylineDouble( spA ) MATRIX *spA; #endif /* End case not STDC */ { int ii, ij, ik, im, ib; int iNoBlocks; int iFirstColumn, iLastColumn; /* [a] : Compute no of blocks of rows to be printed */ /* Item [c] */ if( spA->iNoColumns % ((int) COLUMNS_ACROSS_PAGE) == 0 ) iNoBlocks = (spA->iNoColumns/((int) COLUMNS_ACROSS_PAGE)); else iNoBlocks = (spA->iNoColumns/((int) COLUMNS_ACROSS_PAGE)) + 1; /* [b] : Loop over blocks of rows */ for( ib = 1; ib <= iNoBlocks; ib++ ) { /* Item [d] */ iFirstColumn = (ib-1)*((int) COLUMNS_ACROSS_PAGE) + 1; iLastColumn = MIN( ib*((int) COLUMNS_ACROSS_PAGE) , spA->iNoColumns ); /* [c] : Print title of matrix at top of each block */ /* Item [e] */ if( spA->cpName != (char *) NULL ) printf("\nSkyline MATRIX : \"%s\"\n\n", spA->cpName); else printf("\nSkyline MATRIX : \"UNTITLED\"\n\n"); /* [d] : Label row and column nos */ printf ("row/col "); for( ii = iFirstColumn; ii <= iLastColumn; ii++ ) printf("%3d ", ii); printf("\n"); /* [e] : Print Contents of Matrix */ /* Item [f] */ for( ii = 1; ii <= spA->iNoRows; ii++ ) { printf(" %3d ", ii); for( ij = iFirstColumn; ij <= iLastColumn; ij++) { ik = MIN( ii, ij ); im = MAX( ii, ij ); if((im-ik+1) <= spA->uMatrix.dpp[im-1][0]) printf(" %12.5e ", spA->uMatrix.dpp[ im-1 ][ im-ik+1 ]); else printf(" %12.5e ", 0.0); } printf("\n"); } } } /* * ======================================================================= * Dump Skyline Matrix [iNoRows][iNoColumns] of data type DOUBLE * * Input : MATRIX *spA -- Pointer to matrix data structure. * Output : void * ======================================================================= */ #ifdef __STDC__ void matDumpSkyline( MATRIX *spA ) #else /* Start case not STDC */ void matDumpSkyline( spA ) MATRIX *spA; #endif /* End case not STDC */ { int ii, ij; /* [a] : Check that matrix is defined */ if( spA == NULL ) fatalError("Pointer to matrix is NULL in SkyMatrixDump()", (char *) NULL); /* [b] : Print matrix name */ if( spA->cpName != NULL) printf("\nSkyline Matrix: \"%s\", (%3d X%3d)\n", spA->cpName, spA->iNoRows, spA->iNoColumns); else printf("\nSkyline Matrix: NoName, (%3d X%3d)\n", spA->iNoRows, spA->iNoColumns); /* [c] : Dump contents of skyline matrix */ printf("\nMatrix Array: \n"); for(ii = 1; ii <= spA->iNoColumns; ii++) { printf("Column [%2d], Height = %2d : ", ii, (int) (spA->uMatrix.dpp[ii-1][0])); for (ij = 1; ij <= spA->uMatrix.dpp[ii-1][0]; ij++) printf(" %12.5e", spA->uMatrix.dpp[ii-1][ij]); printf("\n"); } } /* * ================================================================ * Matrix Operations : Addition, Subtraction, Multiplication, Copy. * ================================================================ */ /* * ======================================================================== * matReallocSkyline() : Reallocate Memory for Skyline Matrix. Eliminate * matrix elements that are smaller than MINREAL * * Input : MATRIX *spA -- Pointer to matrix data structure. * Output : MATRIX *spB -- Pointer to new skyline matrix. * ======================================================================== */ #ifdef __STDC__ MATRIX *matReallocSkyline ( MATRIX *spA ) #else /* Start case not STDC */ MATRIX *matReallocSkyline ( spA ) MATRIX *spA; #endif /* End case not STDC */ { MATRIX *spB; int ii, ij, *ld; if( spA == NULL ) fatalError("In matReallocSkyline() : spA = NULL", (char *) NULL); ld = vecAllocInteger( spA->iNoColumns ); for(ii=0; ii < spA->iNoColumns; ii = ii + 1) { for (ij = spA->uMatrix.dpp[ii][0]; ij > 1, abs(spA->uMatrix.dpp[ii][ij]) < MINREAL; ij = ij - 1); ld[ii] = ij; } spB = matAllocSkyline( spA->cpName, DoubleArray, spA->iNoRows, spA->iNoColumns, ld); for(ii = 0; ii < spA->iNoColumns; ii++) for(ij = 1; ij <= spB->uMatrix.dpp[ii][0]; ij++) spB->uMatrix.dpp[ii][ij] = spA->uMatrix.dpp[ii][ij]; return (spB); } /* * ======================================================================= * matAddSkyline() : Add Skyline Matrices * * Input : MATRIX *spA -- Pointer to matrix data structure. * : MATRIX *spB -- Pointer to matrix data structure. * Output : MATRIX *spC -- Pointer to matrix sum * ======================================================================= */ #ifdef __STDC__ MATRIX *matAddSkyline( MATRIX *spA , MATRIX *spB ) #else /* Start case not STDC */ MATRIX *matAddSkyline( spA, spB ) MATRIX *spA, *spB; #endif /* End case not STDC */ { MATRIX *spC; int *ld, min_ldi; int ii, ij; /* [a] : Check for compatible matrix dimensions and data types */ if((spA == NULL) || (spB == NULL)) fatalError("In matAddSkyline() : Either spA = NULL or spB = NULL", (char *) NULL); if(spA->Rep != spA->Rep) fatalError("In matAddSkyline() : spA->Rep != spB->Rep", (char *) NULL ); if(spA->Type != spB->Type) fatalError("In matAddSkyline() : spA->Type != spB->Type", (char *) NULL); /* [b] : Compute Skyline Profile for [A] + [B] */ ld = vecAllocInteger( spA->iNoColumns ); for(ii = 0; ii < spA->iNoColumns; ii++) ld[ii] = (int) MAX (spA->uMatrix.dpp[ii][0], spB->uMatrix.dpp[ii][0]); /* [c] : Allocate Skyline Matrix for [C] : [C] = [A] + [B] */ spC = matAllocSkyline((char *) NULL, DoubleArray, spA->iNoRows, spA->iNoColumns, ld); /* [d] : Compute [C] = [A] + [B] */ for(ii = 0; ii < spA->iNoColumns; ii++) { min_ldi = (int) MIN ( spA->uMatrix.dpp[ii][0], spB->uMatrix.dpp[ii][0]); for(ij = 1; ij <= min_ldi; ij++) spC->uMatrix.dpp[ii][ij] = spA->uMatrix.dpp[ii][ij] + spB->uMatrix.dpp[ii][ij]; if( spA->uMatrix.dpp[ii][0] > min_ldi) for(ij = min_ldi+1; ij <= ld[ii]; ij++) spC->uMatrix.dpp[ii][ij] = spA->uMatrix.dpp[ii][ij]; if( spB->uMatrix.dpp[ii][0] > min_ldi) for(ij = min_ldi+1; ij <= ld[ii]; ij++) spC->uMatrix.dpp[ii][ij] = spB->uMatrix.dpp[ii][ij]; } spC = matReallocSkyline ( spC ); return ( spC ); } /* * ======================================================================= * matSubSkyline() : Subtract two Skyline Matrices * * Input : MATRIX *spA -- Pointer to matrix data structure. * : MATRIX *spB -- Pointer to matrix data structure. * Output : MATRIX *spC -- Pointer to matrix differnce. * ======================================================================= */ #ifdef __STDC__ MATRIX *matSubSkyline( MATRIX *spA , MATRIX *spB ) #else /* Start case not STDC */ MATRIX *matSubSkyline( spA, spB ) MATRIX *spA, *spB; #endif /* End case not STDC */ { MATRIX *spC; int *ld, min_ldi; int ii, ij; /* [a] : Check for compatible matrix dimensions and data types */ if((spA == NULL) || (spB == NULL)) fatalError("In matSubSkyline() : Either spA = NULL or spB = NULL", (char *) NULL); if(spA->Rep != spA->Rep) fatalError("In matSubSkyline() : spA->Rep != spB->Rep", (char *) NULL ); if(spA->Type != spB->Type) fatalError("In matSubSkyline() : spA->Type != spB->Type", (char *) NULL); /* [b] : Compute Skyline Profile for [A] + [B] */ ld = vecAllocInteger( spA->iNoColumns ); for(ii = 0; ii < spA->iNoColumns; ii++) ld[ii] = MAX (spA->uMatrix.dpp[ii][0], spB->uMatrix.dpp[ii][0]); /* [c] : Allocate Skyline Matrix for [C] : [C] = [A] - [B] */ spC = matAllocSkyline(spA->cpName, DoubleArray, spA->iNoRows, spA->iNoColumns, ld); /* [d] : Compute [C] = [A] + [B] */ for(ii = 0; ii < spA->iNoColumns; ii++) { min_ldi = MIN ( spA->uMatrix.dpp[ii][0], spB->uMatrix.dpp[ii][0]); for(ij = 1; ij <= min_ldi; ij++) spC->uMatrix.dpp[ii][ij] = spA->uMatrix.dpp[ii][ij] - spB->uMatrix.dpp[ii][ij]; if( spA->uMatrix.dpp[ii][0] > min_ldi) for(ij = min_ldi+1; ij <= ld[ii]; ij++) spC->uMatrix.dpp[ii][ij] = spA->uMatrix.dpp[ii][ij]; if( spB->uMatrix.dpp[ii][0] > min_ldi) for(ij = min_ldi+1; ij <= ld[ii]; ij++) spC->uMatrix.dpp[ii][ij] = -spB->uMatrix.dpp[ii][ij]; } spC = matReallocSkyline ( spC ); return ( spC ); } /* * ======================================================================= * matMultSkyline() : Multiply two Skyline Matrices. * * Input : MATRIX *spA -- Pointer to matrix data structure [A]. * : MATRIX *spB -- Pointer to matrix data structure [B]. * Output : MATRIX *spC -- Pointer to matrix product [C] = [A]*[B] * * Note : If [A] and [B] are symmetric, then [A]*[B] is symmetric only iff * [A]*[B] = [B]*[A]. In general, this condition will not hold -- so * we store matrix product [ma]*[mb] in Indirect form. * ======================================================================= */ #ifdef __STDC__ MATRIX *matMultSkyline( MATRIX *spA, MATRIX *spB ) #else /* Start case not STDC */ MATRIX *matMultSkyline( spA, spB ) MATRIX *spA; MATRIX *spB; #endif /* End case not STDC */ { MATRIX *spC; int ii, ij, ik, init; int im, in, kka, kkb; /* [a] : Check for compatible matrix dimensions and data types etc... */ if((spA == NULL) || (spB == NULL)) fatalError("In matMultSkyline() : Either spA = NULL or spB = NULL", (char *) NULL); if(spA->Rep != spB->Rep) fatalError("In matMultSkyline() : spA->Rep != spB->Rep", (char *) NULL); if(spA->Type != spB->Type) fatalError("In matMultSkyline() : spA->Type != spB->Type", (char *) NULL); if(spA->iNoColumns != spB->iNoRows) fatalError("In matMultSkyline() : spA->iNoColumns != spB->iNoRows", (char *) NULL); /* [b] : Allocate Memory for Matrix Product */ spC = matAllocIndirect( spA->cpName, DoubleArray, spA->iNoRows, spB->iNoColumns); /* [c] : Compute elements of matrix product */ for( ii=1; ii <= spA->iNoRows; ii++) { for( ij=1; ij <= spB->iNoColumns; ij++) { init = MAX((ii- spA->uMatrix.dpp[ii-1][0]+1),(ij - spB->uMatrix.dpp[ij-1][0]+1)); for(ik = init; ik <= spA->iNoColumns; ik++) { im = MAX(ik,ii); kka = MIN(ik,ii); in = MAX(ik,ij); kkb = MIN(ik,ij); if((im-kka+1) <= spA->uMatrix.dpp[ im ][0] && (in-kkb+1) <= spB->uMatrix.dpp[ in ][0] ) spC->uMatrix.dpp[ ii-1 ][ ij-1 ] += spA->uMatrix.dpp[ im ][ im-kka+1 ]*spB->uMatrix.dpp[ in ][ in-kkb+1 ]; } } } return (spC); } /* * =============================================================================== * matMultSkylineIndirectDouble() : Compute Product of Skyline-Indirect Matrix. * * Input : MATRIX *spA -- Pointer to matrix data structure [A]. * : MATRIX *spB -- Pointer to matrix data structure [B]. * Output : MATRIX *spC -- Pointer to matrix product [C] = [A]*[B] * * Note : If [A] and [B] are symmetric, then [A]*[B] is symmetric only iff * [A]*[B] = [B]*[A]. In general, this condition will not hold -- so * we store matrix product [ma]*[mb] in Indirect form. * =============================================================================== */ #ifdef __STDC__ MATRIX *matMultSkylineIndirectDouble( MATRIX *spA, MATRIX *spB ) #else /* Start case not STDC */ MATRIX *matMultSkylineIndirectDouble( spA, spB ) MATRIX *spA; MATRIX *spB; #endif /* End case not STDC */ { MATRIX *spC; double dAvalue; int ii, ij, ik, init; int im, in; int iMin, iMax; /* [a] : Check for compatible matrix dimensions and data types etc... */ if((spA == NULL) || (spB == NULL)) fatalError("In matMultSkylineIndirect() : Either spA = NULL or spB = NULL", (char *) NULL); if(spA->iNoColumns != spB->iNoRows) fatalError("In matMultSkylineIndirect() : spA->iNoColumns != spB->iNoRows", (char *) NULL); /* [b] : Allocate Memory for Matrix Product */ spC = matAllocIndirect((char *) NULL, DoubleArray, spA->iNoRows, spB->iNoColumns); /* [c] : Compute elements of matrix product */ for( ii=1; ii <= spA->iNoRows; ii++) { for( ij=1; ij <= spB->iNoColumns; ij++) { init = ii - spA->uMatrix.dpp[ii-1][0] + 1; for(ik = init; ik <= spA->iNoColumns; ik++) { im = MIN( ii, ik ); in = MAX( ii, ik ); if((in-im+1) <= spA->uMatrix.dpp[in-1][0]) dAvalue = spA->uMatrix.dpp[ in-1 ][ in-im+1 ]; else dAvalue = 0.0; spC->uMatrix.dpp[ ii-1 ][ ij-1 ] += dAvalue * spB->uMatrix.dpp[ ik - 1 ][ ij - 1 ]; } } } return (spC); } #ifdef __STDC__ MATRIX *matMultIndirectSkylineDouble( MATRIX *spA, MATRIX *spB ) #else /* Start case not STDC */ MATRIX *matMultIndirectSkylineDouble( spA, spB ) MATRIX *spA; MATRIX *spB; #endif /* End case not STDC */ { MATRIX *spC; double dBvalue; int ii, ij, ik, init; int im, in; int iMin, iMax; /* [a] : Check for compatible matrix dimensions and data types etc... */ if((spA == NULL) || (spB == NULL)) fatalError("In matMultSkylineIndirect() : Either spA = NULL or spB = NULL", (char *) NULL); if(spA->iNoColumns != spB->iNoRows) fatalError("In matMultSkylineIndirect() : spA->iNoColumns != spB->iNoRows", (char *) NULL); /* [b] : Allocate Memory for Matrix Product */ spC = matAllocIndirect((char *) NULL, DoubleArray, spA->iNoRows, spB->iNoColumns); /* [c] : Compute elements of matrix product */ for( ii=1; ii <= spA->iNoRows; ii++) { for( ij=1; ij <= spB->iNoColumns; ij++) { init = ii - spB->uMatrix.dpp[ii-1][0] + 1; for(ik = init; ik <= spB->iNoRows; ik++) { im = MIN( ik, ij ); in = MAX( ik, ij ); if((in-im+1) <= spB->uMatrix.dpp[in-1][0]) dBvalue = spB->uMatrix.dpp[ in-1 ][ in-im+1 ]; else dBvalue = 0.0; spC->uMatrix.dpp[ ii-1 ][ ij-1 ] += spA->uMatrix.dpp[ ii - 1 ][ ik - 1 ]* dBvalue; } } } return (spC); } /* * ======================================================================= * matCopySkyline() : Add Skyline Matrices * * Input : MATRIX *spA -- Pointer to matrix data structure. * Output : MATRIX *spCopy -- Pointer to matrix copy * ======================================================================= */ #ifdef __STDC__ MATRIX *matCopySkyline( MATRIX *spA ) #else /* Start case not STDC */ MATRIX *matCopySkyline( spA ) MATRIX *spA; #endif /* End case not STDC */ { MATRIX *spCopy; int *ld, ii, ij; /* [a] : Check input for [ma] */ if(spA == NULL) fatalError("In matCopySkyline() : [A] = NULL", (char *) NULL); /* [b] : Copy Matrix [spCopy] = [spA] */ ld = vecAllocInteger ( spA->iNoColumns); for(ii = 1; ii <= spA->iNoColumns; ii++) ld[ii-1] = spA->uMatrix.dpp[ii-1][0]; spCopy = matAllocSkyline( spA->cpName, DoubleArray, spA->iNoRows, spA->iNoColumns, ld); for(ii = 1; ii <= spA->iNoColumns; ii++) for(ij = 1; ij <= ld[ii-1]; ij++) spCopy->uMatrix.dpp[ii-1][ij] = spA->uMatrix.dpp[ii-1][ij]; 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 *matIndirectToSkyline( MATRIX *spA ) #else /* Start case not STDC */ MATRIX *matIndirectToSkyline( spA ) MATRIX *spA; #endif /* End case not STDC */ { MATRIX *spB; int *ld; int ii, ij; /* [a] : Check input matrix [A] is defined and symmetric */ if(spA == NULL) fatalError("In matIndirectToSkyline() : spA == NULL", (char *) NULL); for(ii = 1; ii <= spA->iNoColumns-1; ii++) { for(ij = 0; ij <= ii-1; ij++) if(abs( spA->uMatrix.dpp[ij][ii] - spA->uMatrix.dpp[ii][ij]) > MINREAL) fatalError("In matIndirectToSkyline() : [A] not symmetric", (char *) NULL); } /* [c] : Compute required Ld[] matrix */ ld = vecAllocInteger( spA->iNoColumns ); for(ii = 1; ii <= spA->iNoColumns; ii++) { for(ij = 1; ij < ii, ABS(spA->uMatrix.dpp[ij-1][ii-1]) < MINREAL; ij++); ld[ii-1] = ii-ij+1; } spB = matAllocSkyline(spA->cpName, DoubleArray, spA->iNoRows, spA->iNoColumns, ld); switch ((int) spA->Type ) { case DoubleArray: for(ii = 1; ii <= spA->iNoRows; ii++) for(ij = 1; ij <= ld[ii-1]; ij++) spB->uMatrix.dpp[ii-1][ij] = spA->uMatrix.dpp[ii-ij][ii-1]; break; default: fatalError("In matIndirectToSkyline() : Undefined matrix type", (char *) NULL); } return ( spB ); } /* * =================================================================== * 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 */ { MATRIX *spB = NULL; return (spB); } /* * ======================================================================= * matScaleSkyline() : 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 *matScaleSkyline( MATRIX *spA, double dScale) #else /* Start case not STDC */ MATRIX *matScaleSkyline( 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 matIndirectToSkyline() : [A] == NULL", (char *) NULL); spB = matCopySkyline( spA ); for(ii = 1; ii <= spA->iNoColumns; ii++) for(ij = 1; ij <= spA->uMatrix.dpp[ii-1][0]; ij++) spB->uMatrix.dpp[ii-1][ij] = dScale * spA->uMatrix.dpp[ii-1][ij]; return ( spB); } /* * ========================================================================= * matLUDecompositionSkyline() : [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 *matLUDecompositionSkyline( MATRIX *spA ) #else /* Start case not STDC */ MATRIX *matLUDecompositionSkyline( spA ) MATRIX *spA; #endif /* End case not STDC */ { 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 */ switch((int) spA->Type ) { case DoubleArray: for(ii = 0; ii < spA->iNoColumns; ii++) { for(ij = spA->uMatrix.dpp[ii][0]; ij > 1; ij = ij - 1) { min_col = MIN((spA->uMatrix.dpp[ii][0]), (spA->uMatrix.dpp[ii-ij+1][0]+ij-1)); for(ik=ij+1; ik <= min_col; ik++) spA->uMatrix.dpp[ii][ij] -= spA->uMatrix.dpp[ii][ik] * spA->uMatrix.dpp[ii-ij+1][ik-ij+1] * spA->uMatrix.dpp[ii-ik+1][1]; spA->uMatrix.dpp[ii][ij] /= spA->uMatrix.dpp[ii-ij+1][1]; } for(ij = spA->uMatrix.dpp[ii][0]; ij > 1; ij--) spA->uMatrix.dpp[ii][1] -= spA->uMatrix.dpp[ii][ij] * spA->uMatrix.dpp[ii][ij] * spA->uMatrix.dpp[ii-ij+1][1]; if(abs(spA->uMatrix.dpp[ii][1]) <= MINREAL) fatalError("In LUDecompositionSkyline() : Singular Matrix ", (char *) NULL); } break; default: fatalError("In LUDecompositionSkyline() : Undefined Matrix Data Type", (char *) NULL); } return ( spA ); } /* * ===================================================================================== * matLUBacksubstitutionSkyline() : 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 *matLUBacksubstitutionSkyline( MATRIX *spA, MATRIX *spB) #else /* Start case not STDC */ MATRIX *matLUBacksubstitutionSkyline( spA, spB ) MATRIX *spA; MATRIX *spB; #endif /* End case not STDC */ { int ii, ij, ik; #ifdef MYDEBUG printf("*** Enter matLUBacksubstitutionSkyline()\n"); #endif /* [a] : Check compatibility of input matrices */ if(spA == NULL) { fatalError("In matLUBacksubstitutionSkyline() : spA == NULL", (char *) NULL); } if(spA->iNoColumns != spB->iNoRows) { fatalError("In matLUBacksubstitutionSkyline() : Dimensions of", "A[][] and b[] are incompatible ", (char *) NULL); } /* [b] : Compute forward substitution */ for(ij=0; ij < spB->iNoColumns; ij++) for(ii = 0; ii < spB->iNoRows; ii++) { for(ik = spA->uMatrix.dpp[ii][0]; ik > 1; ik--) spB->uMatrix.dpp[ii][ij] -= spB->uMatrix.dpp[ii-ik+1][ij] * spA->uMatrix.dpp[ii][ik]; } for(ij = 0; ij < spB->iNoColumns; ij++) for(ii=0; ii < spB->iNoRows; ii++) spB->uMatrix.dpp[ii][ij] /= spA->uMatrix.dpp[ii][1]; /* [c] : Compute backward substitution */ for(ij = 0; ij < spB->iNoColumns; ij++) for(ii = spB->iNoRows-1; ii >= 0; ii--) for(ik = ii+1; ik <= spB->iNoRows-1; ik++) if((ik+1-ii) <= spA->uMatrix.dpp[ik][0] ) spB->uMatrix.dpp[ii][ij] -= spB->uMatrix.dpp[ik][ij] * spA->uMatrix.dpp[ik][ik+1-ii]; return ( spB ); } /* * =================================================================== * matInverseSkyline() : 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 *matInverseSkyline ( MATRIX *spA ) #else /* Start case not STDC */ MATRIX *matInverseSkyline ( spA ) MATRIX *spA; #endif /* End case not STDC */ { MATRIX *p; MATRIX *pn; int ii, ij; /* [a] : Check Input */ if( spA == NULL ) fatalError("In matInverseSkyline () : [ma] == NULL \n", (char *) NULL); /* [b] : Compute LU Decomposition and column vectors of inverse */ p = matCopySkyline( spA ); pn = matAllocIndirect( spA->cpName, Indirect, spA->iNoRows, spA->iNoColumns); for(ii = 0; ii < spA->iNoRows; ii++) for(ij = 0; ij < spA->iNoColumns; ij++) if(ii == ij) pn->uMatrix.dpp[ii][ij] = 1.0; else pn->uMatrix.dpp[ii][ij] = 0.0; p = matLUDecompositionSkyline(p); pn = matLUBacksubstitutionSkyline( p, pn); matFreeSkyline(p); p = matIndirectToSkyline(pn); matFree (pn); p = matReallocSkyline(p); return (p); } /* * ==================================================================== * matSubspace() : 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 : Mark Austin November 1993 * ==================================================================== */ enum { MaxSubspaceIterations = 50 }; #define TOLERANCE 0.00000001 #ifdef __STDC__ void matSubspace( MATRIX *spA, MATRIX *spB, MATRIX *spEigenvalue, MATRIX *spEigenvector, int iNoEigen ) #else /* Start case not STDC */ void matSubspace( spA, spB, spEigenvalue, spEigenvector , iNoEigen) MATRIX *spA, *spB; MATRIX *spEigenvalue; MATRIX *spEigenvector; int iNoEigen; #endif /* End case not STDC */ { MATRIX *spEigenvectorQ, *spEigenvalueOld; MATRIX *spAwork, *spBwork; MATRIX *spTemp1, *spA_LU; MATRIX *spY, *spYnew, *spV, *spVk; double dEigenvalueOld, dConvergence; double dEigenvalue, dMaxValue; int iSize, ii, ij, ik; int iMin, iMax, iLoop; int iSubspaceConvergence; /* [a] : Check Input and Allocate Working Matrices */ if( spA == NULL || spB == NULL ) fatalError("In matSubspace() : Pointer spA == NULL or spB == NULL", (char *) NULL); if( spA->Type != DoubleArray || spB->Type != DoubleArray ) fatalError("In matSubspace() : spA->Type (or spB->Type) != DoubleArray", (char *) NULL); if( spA->Rep != Skyline || spB->Rep != Skyline ) fatalError("In matSubspace() : spA->Rep != Skyline or spB->Rep != Skyline", (char *) NULL); /* [b] : Initialize Starting Eigenvectors [X_{k}] */ iSize = spA->iNoRows; for(ii = 1; ii <= iNoEigen; ii = ii + 1) spEigenvector->uMatrix.dpp[ii-1][ii-1] = 1.0; /* [c] : Loops for Subspace Iteration */ spAwork = matAllocIndirect("[A_{k+1}]", DoubleArray, iNoEigen, iNoEigen); spBwork = matAllocIndirect("[B_{k+1}]", DoubleArray, iNoEigen, iNoEigen); spEigenvectorQ = matAllocIndirect("[Q]", DoubleArray, iNoEigen, iNoEigen); spVk = matAllocIndirect("[V_{k}]", DoubleArray, iSize, iNoEigen); spY = matAllocIndirect("[Y_{k}]", DoubleArray, iSize, iNoEigen); spYnew = matAllocIndirect("[Y_{k+1}]", DoubleArray, iSize, iNoEigen); spEigenvalueOld = matCopy( spEigenvalue ); iLoop = 0; iSubspaceConvergence = False; while( iLoop <= (int) MaxSubspaceIterations && iSubspaceConvergence == False ) { iLoop = iLoop + 1; /* [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.dpp[ 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.dpp[iMax-1][0]) { spY->uMatrix.dpp[ ii-1 ][ ij-1 ] += spB->uMatrix.dpp[ iMax-1 ][ iMax-iMin+1 ] * spEigenvector->uMatrix.dpp[ ik-1 ][ ij-1 ] ; } } } /* [c.2] : Solve [A][V_{k+1}] = [Y_{k}] */ if(iLoop == 1) { spA_LU = matLUDecompositionSkyline(matCopy(spA)); } for(ii = 1; ii <= iSize; ii = ii + 1) for(ij = 1; ij <= iNoEigen; ij = ij + 1) spVk->uMatrix.dpp[ ii-1 ][ ij-1 ] = spY->uMatrix.dpp[ ii-1 ][ ij-1 ]; if(iLoop != 1) { matFree (spTemp1); } spV = matLUBacksubstitutionSkyline( 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.dpp[ 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.dpp[iMax-1][0]) { spYnew->uMatrix.dpp[ ii-1 ][ ij-1 ] += spA->uMatrix.dpp[ iMax-1 ][ iMax-iMin+1 ] * spV->uMatrix.dpp[ 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.dpp[ ii-1 ][ ij-1 ] = 0.0; for(ik = 1; ik <= iSize; ik = ik + 1) spAwork->uMatrix.dpp[ ii-1 ][ ij-1 ] += spV->uMatrix.dpp[ ik-1 ][ ii-1 ] * spYnew->uMatrix.dpp[ 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.dpp[ 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.dpp[iMax-1][0]) { spYnew->uMatrix.dpp[ ii-1 ][ ij-1 ] += spB->uMatrix.dpp[ iMax-1 ][ iMax-iMin+1 ] * spV->uMatrix.dpp[ 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.dpp[ ii-1 ][ ij-1 ] = 0.0; for(ik = 1; ik <= iSize; ik = ik + 1) spBwork->uMatrix.dpp[ ii-1 ][ ij-1 ] += spV->uMatrix.dpp[ ik-1 ][ ii-1 ] * spYnew->uMatrix.dpp[ ik-1 ][ ij-1 ]; } /* [c.7] : Solve [A_{k+1}].[Q_{k+1}] = [B_{k+1}].[Q_{k+1}].[Lambda_{k+1}] */ matGeneralisedEigen( spAwork, spBwork, spEigenvalue, spEigenvectorQ ); /* [c.8] : Check Convergence Criteria */ iSubspaceConvergence = True; for( ii = 1; ii <= spEigenvalue->iNoRows; ii = ii + 1) { dEigenvalue = spEigenvalue->uMatrix.dpp[ii-1][0]; dEigenvalueOld = spEigenvalueOld->uMatrix.dpp[ii-1][0]; dConvergence = fabs((dEigenvalueOld - dEigenvalue))/dEigenvalue; if(dConvergence > TOLERANCE) { iSubspaceConvergence = False; } } matFree ( spEigenvalueOld ); spEigenvalueOld = matCopy( spEigenvalue ); /* [c.9] : Update [X_{k+1}] = [V_{k+1}].[Q_{k+1}] */ spTemp1 = matMult( 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.dpp[ ii-1 ][ ij-1 ]) >= ABS(dMaxValue)) dMaxValue = spTemp1->uMatrix.dpp[ ii-1 ][ ij-1 ]; for(ii = 1; ii <= iSize; ii = ii + 1) spEigenvector->uMatrix.dpp[ ii-1 ][ ij-1 ] = spTemp1->uMatrix.dpp[ ii-1 ][ ij-1 ]/dMaxValue; } } /* [d] : Summarize Algorithm Performance : clean up */ printf("\n*** Subspace Iteration converged in %2d iterations \n", iLoop); matFree ( spTemp1 ); matFree ( spAwork ); matFree ( spBwork ); matFree ( spEigenvectorQ ); matFree ( spEigenvalueOld ); matFree ( spYnew ); matFree ( spY ); }