matrix.c

Embedded magazine source code. this is the source code in the year 1989

    } /* while( pChain != NULL */

} /* remove_chain */


/*** transpose() ***********************************************************
 *
 *   Mike Foley
 *   January 21, 1989
 *
 *   FUNCTIONAL SUMMARY :
 *	Perform an in place transposition of the given matrix.
 *
 *   INPUT :
 *	ret_condition = transpose( &pMatrix )
 *
 *	    pMatrix : Pointer to the matrix to be tranposed.
 *
 *   OUTPUT :
 *	0 : Function performed as expected.
 *	1 : Error, out of memory.
 *
 *   OPERATIONS :
 *	From the matrix type determine the number of chains, and the maximum
 *	    length of each chain.
 *	Dimension the resulting matrix to the proper dimension for the
 *	    transposed matrix.	This means swap the row and column dimensions
 *	    of the given matrix.
 *	For each position in the original matrix, save it in the new matrix
 *	    in the reversed position ( column, row ).
 *	Free the storage used by the original matrix.
 *	Point the old matrix pointer to the resultant matrix.
 *
 ***************************************************************************/
int transpose( MATRIX **mat )
{
    MATRIX *pResult;
    MATRIX *pMat;
    BYTE bLinks;
    BYTE bLength;
    BYTE bType;
    BYTE bCurElem;
    BYTE bTransRow;
    BYTE bTransCol;
    BYTE bRow;
    BYTE bCol;
    double dVal;

    pMat = *mat;
    bType = pMat->bType & LOW;
    bRow = pMat->bRow;
    bCol = pMat->bCol;

    if( bType == ROW ) {
	bLinks = pMat->bRow;
	bLength = pMat->bCol;
	bTransRow = bLength;
	bTransCol = bLinks;
    } else {
	bLinks = pMat->bCol;
	bLength = pMat->bRow;
	bTransRow = bLinks;
	bTransCol = bLength;
    }

    if( ( pResult = dimension( bTransRow, bTransCol, bType ) ) == NULL )
	return( 1 );

    while( bRow-- > 0 ) {

	bCurElem = bCol;
	while( bCurElem-- > 0 ) {

	    out_matrix( pMat, bRow, bCurElem, &dVal );
	    in_matrix( pResult, bCurElem, bRow, dVal );

	} /* while( bCurElem-- > 0 ) */

    } /* while( bRow-- > 0 ) */

    remove_matrix( pMat );
    *mat = pResult;
    return( 0 );

} /* transpose */


/*** mat_print() ********************************************************
 *
 *   Mike Foley
 *   January 21, 1989
 *
 *   FUNCTIONAL SUMMARY :
 *	Display onscreen, the entire matrix.  This includes zero elements,
 *    not just elements actually stored in the matrix.
 *
 *   INPUT :
 *	ret_condition = mat_print( pMatrix )
 *
 *	    pMatrix : Pointer to matrix to be displayed onscreen.
 *
 *   OUTPUT :
 *	0 : Matrix was displayed as expected.
 *	1 : Error, given matrix was NULL
 *
 *   OPERATIONS :
 *	Ensure matrix isn't NULL
 *	Use the matrices dimensions to loop through each coordinate.
 *	Display whatever is in that position in the matrix.
 *	Only print newlines after each each row is complete.
 *	NOTE - Formattint isn't performed.  If a matrix has many columns,
 *	    the rows will flow onto the next screen line.
 *
 ***************************************************************************/
int mat_print( MATRIX *pMatrix )
{
    BYTE bRow;
    BYTE bCol;
    BYTE i;
    BYTE j;
    double dTemp;

    if( pMatrix == NULL )
	return( 1 );

    bRow = pMatrix->bRow;
    bCol = pMatrix->bCol;

    for( i = 0; i < bRow; i++ ) {

	for( j = 0; j < bCol; j++ ) {

	    out_matrix( pMatrix, i, j, &dTemp );
	    printf( "%f ", dTemp );

	} /* for( j = 0; j < bCol; j++ ) */

	printf( "\n" );

    } /* for( i = 0; i < bRow; i++ ) */

    return( 0 );

} /* mat_print */



/*** solve() ********************************************************
 *
 *   Mike Foley
 *   January 21, 1989
 *
 *   FUNCTIONAL SUMMARY :
 *	Solve the equation Ax = b for x; given the matrices A and b.
 *
 *   INPUT :
 *	x = solve( pA, pb )
 *
 *	    A : Pointer to n x n matrix.
 *	    b : Pointer to n x 1 matrix ( vector ).
 *
 *   OUTPUT :
 *	MATRIX * : Solution of equation Ax = b.
 *	NULL : Error condition.
 *	       Out of memory.
 *	       A or b matrix is NULL.
 *	       A matrix not square.
 *	       Dimensions of A and b not compatible.
 *
 *   OPERATIONS :
 *	Determine valididy of matrices.
 *	The A matrix must be in LU format.  If it hasn't already been
 *	    decomposed, have it done now.  Save this representation in a
 *	    local variable.  If the matrix is already decomposed, point the
 *	    local variable to this matrix, and set the NOREMOVE flag in it's
 *	    type field.
 *	Perform forward and backward substitution using the LU factorization
 *	    and the b matrix.  These steps have been placed in separate
 *	    functions.	If it is know that a matrix is triangular, the
 *	    proper function can be called directly, eliminating the need to
 *	    use solve.	This will speed solutions for this special case of
 *	    A matrices.
 *	Remove storage used for tempory matrices.  If the matrix A was factored
 *	    for this routine, remove it's storage.
 *	Return the solution x.
 *
 ***************************************************************************/
MATRIX *solve( MATRIX *pA, MATRIX *pb )
{
    MATRIX *pMiddle;
    MATRIX *pX;
    MATRIX *pLU;

    if( ( pA == NULL ) || ( pb == NULL ) )
	return( NULL );
    if( ( pA->bCol != pb->bRow ) || ( pA->bCol != pA->bRow ) )
	return( NULL );
    if( pb->bCol != 1 )
	return( NULL );

    if( pA->bType & LU ) {
	pLU = pA;
	pLU->bType += NOREMOVE;
    } else
	pLU = decompose( pA );

    pMiddle = forward_sub( pLU, pb );
    pX = back_sub( pLU, pMiddle );

    remove_matrix( pMiddle );
    if( pLU->bType & NOREMOVE )
	pLU->bType -= NOREMOVE;
    else
	remove_matrix( pLU );

    return( pX );

} /* solve */


/*** forward_sub() ********************************************************
 *
 *   Mike Foley
 *   January 21, 1989
 *
 *   FUNCTIONAL SUMMARY :
 *	Perform forward substitution on the given matrices.
 *
 *   INPUT :
 *	result_matrix = forward_sub( pMatrix, pVector )
 *
 *	    pMatrix : Pointer to n x n LU factored matrix.
 *	    PVector : Pointer to n x 1 constant Vector.
 *
 *   OUTPUT :
 *	MATRIX * : Solution of forward substitution n x 1 vector.
 *
 *   OPERATIONS :
 *	Insure matrices are dimensioned properly for funcion.
 *	Dimension resultant matrix to same dimensions and type as the input
 *	    vector.
 *	Perform the forward substitution operation.
 *	Return resultant vector.
 *
 ***************************************************************************/
MATRIX *forward_sub( MATRIX *mat, MATRIX *vector )
{
    MATRIX *pResult;
    MATRIX_ELEM *pNext;
    BYTE bVRow;
    BYTE bVCol;
    BYTE bMax;
    double dSum;
    double dTemp1;
    double dTemp2;
    USHORT i;

    bVCol = vector->bCol;
    bVRow = vector->bRow;
    if( bVCol != 1 )
	return( NULL );
    if( !(vector->bType & COL ) )
	change_rep( &vector );

    if( ( pResult = dimension( bVRow, bVCol, COL ) ) == NULL )
	return( NULL );

    pNext = ( *vector->pFirst )[ 0 ];
    while( pNext != NULL ) {

	dSum = 0.0;
	bMax = pNext->index;
	for( i = 0; i < bMax; i++ ) {

	    out_matrix( mat, bMax, i, &dTemp1 );
	    out_matrix( vector, i, 0, &dTemp2 );
	    dSum += dTemp1 * dTemp2;

	} /* for( i = 0; i < bMax; i++ ) */
	dSum = pNext->elem - dSum;
	in_matrix( pResult, bMax, 0, dSum );
	pNext = pNext->pNext;

    } /* while( pNext != NULL ) */

    return( pResult );

} /* forward_sub */


/*** back_sub() ********************************************************
 *
 *   Mike Foley
 *   January 21, 1989
 *
 *   FUNCTIONAL SUMMARY :
 *	Perform backward substitution on given matrices.
 *
 *   INPUT :
 *	result_matrix = back_sub( pMatrix, pVector )
 *
 *	    pMatrix : Pointer to n x n LU factored matrix.
 *	    pVector : Pointer to n x 1 constant matrix.
 *
 *   OUTPUT :
 *	MATRIX * : Result of backward substitution n x 1 COL matrix.
 *	NULL : Error condition :
 *	       Out of memory
 *	       Matrices dimensioned improperly.
 *	       One, or both matrices NULL.
 *
 *   OPERATIONS :
 *	Ensure matrices are of proper dimension and representations.
 *	Dimension resultant matrix of dimension equal to input vector.
 *	Perform backward substitution operation.
 *	Return resultant matrix.
 *
 ***************************************************************************/
MATRIX *back_sub( MATRIX *mat, MATRIX *vector )
{
    MATRIX *pResult;
    BYTE bVRow;
    BYTE bVCol;
    BYTE bMax;
    double dSum;
    double dTemp1;
    double dTemp2;
    USHORT i;

    bVCol = vector->bCol;
    bVRow = vector->bRow;
    if( bVCol != 1 )
	return( NULL );
    if( !(vector->bType & COL ) )
	change_rep( &vector );

    if( ( pResult = dimension( bVRow, bVCol, COL ) ) == NULL )
	return( NULL );

    bMax = bVRow;
    while( bMax-- > 0 ) {

	dSum = 0.0;
	i = bVRow;
	while( --i > bMax ) {
	    out_matrix( mat, bMax, i, &dTemp1 );
	    out_matrix( pResult, i, 0, &dTemp2 );
	    dSum += dTemp1 * dTemp2;
	}
	out_matrix( vector, bMax, 0, &dTemp1 );
	dTemp1 -= dSum;
	out_matrix( mat, bMax, bMax, &dTemp2 );
	dTemp2 = dTemp1 / dTemp2;
	in_matrix( pResult, bMax, 0, dTemp2 );
    }

    return( pResult );

} /* back_sub */


/*** decompose() ********************************************************
 *
 *   Mike Foley
 *   January 21, 1989
 *
 *   FUNCTIONAL SUMMARY :
 *	Determine the LU factorization of the given matrix.  Do not
 *   overwrite the input matrix.
 *
 *   INPUT :
 *	LU_factorization = decompose( pMatrix )
 *
 *	    pMatrix : Pointer to n x n input matrix.
 *
 *   OUTPUT :
 *	MATRIX * : Pointer to n x n matrix containing LU factors of input matrix
 *	NULL : Error condition :
 *	       Input matrix NULL.
 *	       Input matrix not square.
 *	       Out of memory.
 *
 *   OPERATIONS :
 *	Ensure input matrix isn't NULL.
 *	Copy input matrix to output matrix.  This copied matrix will become
 *	    the LU factorization.
 *	Only perform LU factorization, if the input matrix isn't already
 *	    factored.
 *	Make sure matrix is square.
 *	If resultant matrix isn't already in COL representation, change it.
 *	Perform a columnwise LU factorization on the resultant matrix.	The
 *	    input matrix, which is in ROW representation, is used to find
 *	    where elements are.
 *	    The factor which will be stored in the lower portion ( L ) of the
 *		resultant matrix is calculated first.  This is the diagonal
 *		element with indices equal to the current column being factored,
 *		divided by the current element.
 *	    The rest of the row which the current element is in, must be
 *		multiplied by the factor, and this product is subtracted from
 *		the element in the same column, but in the row of the
 *		diagonal element used to form the factor.
 *	    These elements are stored in the resultant matrix as they are formed.
 *	    After the row has been reduced, the factor is stored in the L
 *		position, which is where it originally started.
 *	Tag the new matrix as being a LU factorization.
 *	Return LU factorization of the input matrix.
 *
 ***************************************************************************/
MATRIX *decompose( MATRIX *pMatrix )
{
    MATRIX *pResult;
    MATRIX_ELEM *pNext;
    MATRIX_ELEM *pRow;
    MATRIX_ELEM *pCurCol;
    double dFactor;
    double dTemp;
    USHORT i;
    BYTE bType;
    BYTE bRow;
    BYTE bCol;

    if( pMatrix == NULL )
	return( NULL );
    bType = pMatrix->bType;

    if( ( pResult = copy_matrix( pMatrix, COL ) ) == NULL )
	    return( NULL );

    if( !( bType & LU ) ) {

	if( pMatrix->bCol != pMatrix->bRow )
	    return( NULL );
	if( ( pMatrix->bType & LOW ) != ROW )
	    change_rep( &pMatrix );

	bRow = pMatrix->bRow;
	for( i = 0; i < bRow; i++ ) {
	    pNext = ( *pResult->pFirst )[ i ];
	    pRow = ( *pMatrix->pFirst )[ i ];
	    while( pNext != NULL ) {

		bCol = pNext->index;
		if( bCol > i ) {

		    pCurCol = find_location( pRow, i );
		    dTemp = pCurCol->elem;
		    dFactor = pNext->elem / dTemp;
		    pCurCol = pCurCol->pNext;
		    while( pCurCol != NULL ) {

			out_matrix( pResult, bCol, pCurCol->index, &dTemp );
			dTemp -= dFactor * pCurCol->elem ;
			in_matrix( pResult, bCol, pCurCol->index, dTemp );
			pCurCol = pCurCol->pNext;

		    } /* while( pCurCol != NULL ) */

		in_matrix( pResult, bCol, i, dFactor );
		} /* if( bCol < i ) */

		pNext = pNext->pNext;

	    } /* while( pNext != NULL ) */

	} /* for( i = 0; i < pRow; i++ ) */

	pResult->bType += LU;

    } /* if( !( bType & LU ) )	*/

    return( pResult );

} /* decompose */


/*** copy_matrix() ********************************************************
 *
 *   Mike Foley
 *   January 21, 1989
 *
 *   FUNCTIONAL SUMMARY :
 *	Copy the given matrix into a new matrix with the given representation.
 *
 *   INPUT :
 *	new_matrix = copy_matrix( pMatrix, bType )
 *
 *	    pMatrix : Pointer to matrix to be copied.
 *	    bType : Type which new matrix is to be.  Only lower four bits
 *		    can be specified.
 *
 *   OUTPUT :
 *	MATRIX * : Copy of given matrix in desired representation.
 *	NULL : Error condition :
 *	       Out of memory.
 *	       Original matrix was NULL, not truely an error, but a copy of
 *		  the input matrix.
 *
 *   OPERATIONS :
 *	Dimension the resultant matrix to same dimension as input matrix, and
 *	    of type passed to the function.
 *	Scan each row of the input matrix.  For each element in the input
 *	    matrix, copy it to the result matrix at the same postition.
 *	Return the resultant matrix.
 *
 ***************************************************************************/
MATRIX *copy_matrix( MATRIX *pMatrix, BYTE bType )
{
    MATRIX *pResult;
    MATRIX_ELEM *pNext;
    BYTE bRow;
    BYTE bCol;
    BYTE bCurCol;
    double dVal;
    USHORT i;

    bRow = pMatrix->bRow;
    bCol = pMatrix->bCol;
    if( (pResult = dimension( bRow, bCol, bType ) ) == NULL )
	return( NULL );

    bType = pMatrix->bType & LOW;
    for( i = 0; i < bRow; i++ ) {

	pNext = ( *pMatrix->pFirst )[ i ];
	while( pNext != NULL ) {

	    bCurCol = pNext->index;
	    dVal = pNext->elem;
	    if( bType == ROW )
		in_matrix( pResult, i, bCurCol, dVal );
	    else
		in_matrix( pResult, bCurCol, i, dVal );
	    pNext = pNext->pNext;

	} /* while( pNext != NULL ) */

    } /* for( i = 0; i < bRow; i++ ) */

    return( pResult );

} /* copy_matrix */