/* qinv.c - Code file for some quick matrix inverse routines.
* © Copyright 1999 by John Halleck
* All Rights Reserved
*/
/* Version of September 8th, 1999 */
#include "errors.h"
/* We use the package wide standard error codes. */
#include "vec.h"
/* We make use of the vector routines. */
#include "mat.h"
/* We make use of the matrix routines. */
/* #include "matdebug.h" */
/* DEBUG */
/* =============== Inverse without an attempt to pivot ====================== */
/* Non-pivoting inverse. Clearly not a general inverse */
static error invnp (int size, matrix result, matrix working) {
/* Result should either be the identity matrix, or should be
* the result of pivoting the identity matrix to match
* whatever pivoting was done to the working array.
*/
error status; /* Place to capture status of routines we call */
int i, j; /* Indices for loops. (Yes, I used to program in ForTran) */
double *diag, dvalue; /* A diagonal element */
double *leading, lvalue; /* leading non-zero value on a row */
double *workrow, *resultrow; /* The row we are working on. */
double *worksub, *resultsub; /* The other rows we are clearing out. */
/* A stable way to invert a positive definite matrix
* is to use elementary matrix transformations to reduce it
* the identity, while using the same transformations to transform
* the identity to the inverse.
*/
/* Oddly enough, this routine should be much cleaner in Pascal */
if (!result || !working) return ERRnil; /* Arrays must exist */
if (size<1) return ERRsize; /* with reasonable sizes */
if (result==working) return ERRsame; /* And can't be the same arrays. */
workrow = (double *) working; resultrow = (double *) result;
diag = workrow;
for (i=0; i<size; i++) { /* clear up each row. */
/* Scale the row by the leading non-zero entry. */
if (*diag == 0.0) return ERRnumeric;
/* This shouldn't be possible for a positive definite matrix. */
else if (*diag != 1.0) {
dvalue = 1.0 / *diag;
if ((status = vecscl(size, workrow, workrow, dvalue)))
return status;
if ((status = vecscl(size, resultrow, resultrow, dvalue)))
return status;
}
/* Add appropriate scaled version of this row to the other rows
* in order to zero items.
*/
worksub = (double *) working; resultsub = (double *) result;
leading = i + (double *) working;
for (j=0; j<size; j++) {
if (i != j) { /* Don't do it to this row itself. */
if (*leading != 0.0) {
lvalue = - *leading;
if ((status = vecadds (size, worksub, worksub, workrow, lvalue)))
return status;
if ((status = vecadds (size, resultsub, resultsub, resultrow,
lvalue)))
return status;
}
}
leading = leading + size; worksub += size; resultsub += size;
}
workrow += size; resultrow += size; /* drop down a row */
diag += size + 1; /* The diagonal is one row down and one item over */
}
return 0;
}
/* ==================== Inverse without need of pivoting ==================== */
/* Inverse of a positive definite matrix. result and working must be
* distinct arrays.
* Note that positive definite matricies do not require pivoting or
* other rearrangement to form the inverse.
*
* * This is NOT a GENERAL matrix inverse *
*
* This does not make use of the fact that all the positive definite
* matrices in this package are symmetrical
*/
error invpd (int size, matrix result, matrix given, matrix working) {
/* A stable way to invert a positive definite matrix
* is to use elementary matrix transformations to reduce it
* the identity, while using the same transformations to transform
* the identity to the inverse.
* Of course, that would distroy the original, which we probably
* don't want to do...
* So we ask for a matrix for working space.
* Note that setting given and working to the same array means
* we can trash given, and it can run a bit faster.
*/
/* Oddly enough, this routine should be much cleaner in Pascal */
error status; /* Place to capture status of routines we call */
if (!result || !given || !working)
return ERRnil; /* Arrays must exist */
if (size<1)
return ERRsize; /* with reasonable sizes */
if (result==working)
return ERRsame; /* And result and working have to be distinct arrays */
/* Put array into working memory, if necessary */
if (working!=given) if ((status = matcpy (size, size, working, given)))
return status;
/* Get us an identity array to play with. */
if ((status = matidn (size, result))) return status;
/* Return the non-pivoting inverse */
return invnp (size, result, working);
}
/* ==================== Inverse with pivoting =============================== */
/* Invert non-singular array */
error invns (int size, matrix result, matrix given, matrix working) {
/* One can increase the classes of matrix that one can invert
* if one does matrix pivoting.
* This tries the matrix pivot before trying it..
*/
/* Oddly enough, this routine should be much cleaner in Pascal */
/* This routine only does partial pivoting. */
double *workrow, *testvalue, testmax, currentmax;
double *resrow, *rowofmax, *rowofres, *resrowofmax, *rowoftest;
int i, j;
error status;
if (!result || !given || !working)
return ERRnil; /* Arrays must exist */
if (size<1)
return ERRsize; /* with reasonable sizes */
if (result==working)
return ERRsame; /* And can't be the same arrays. */
/* Put array into working memory, if necessary */
if (working!=given) if ((status = matcpy (size, size, working, given)))
return status;
/* Get us an identity array to play with. */
if ((status = matidn (size, result))) return status;
/* Pivot the arrays. */
/* We are doing row swaps to pivot, but we are doing them to
* both arrays.
*/
workrow = (double *) working; resrow = (double *) result;
for (i=0; i<size-1; i++) {
rowofmax = workrow; rowofres = resrow;
/* pointers to current rows in each array */
rowoftest = workrow; /* Row we are testing */
testvalue = workrow + i;
currentmax = *testvalue; if (currentmax < 0) currentmax = -currentmax;
/* The value of the i'th column of the test row is our guess for max */
/* Find the row with the maximum value */
for (j=i+1; j<size; j++) {
rowoftest += size; rowofres += size;
testvalue += size; testmax = *testvalue;
testmax = *testvalue; if (testmax<0) testmax = -testmax;
/* If this row has a bigger value in this column than the prior
* max, it is the new max.
*/
if (testmax > currentmax) {
currentmax = testmax; /* Max value so far */
rowofmax = rowoftest; resrowofmax = rowofres;
/* Row (in each array) where we found the max values */
}
}
/* (The max value had better be greater than 0) */
if (currentmax == 0.0) return ERRnumeric;
/* If the max value is not the current row, swap them. */
if (rowofmax != workrow) { /* swap */
if ((status = vecswp (size, workrow, rowofmax)))
return status;
if ((status = vecswp (size, resrow, resrowofmax)))
return status;
}
workrow += size; resrow += size;
}
/* Return the non-pivoting inverse */
return invnp (size, result, working);
}
This page is http://www.cc.utah.edu/~nahaj/cave/survey/code/c/qinv.c.html
© Copyright 2000 by John Halleck, All Rights Reserved.
This snapshot was last modified on August 23rd, 2000
And the underlying file was last modified on May 11th, 2000