home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
vis-ftp.cs.umass.edu
/
vis-ftp.cs.umass.edu.tar
/
vis-ftp.cs.umass.edu
/
pub
/
Software
/
ASCENDER
/
ascender.tar.Z
/
ascender.tar
/
Triangulate
/
array_utils.c
next >
Wrap
C/C++ Source or Header
|
1995-04-13
|
13KB
|
556 lines
/* @(#)array_utils.c 1.15 12/20/94 */
/*>>
Utility routines for the Numerical recipes cookbook programs
<<*/
#include "cvar.h"
#include "rh_util.h"
#include "alloc.h"
#include <math.h>
#include "newgaussj.s"
#include "newgaussj2.s"
FUNCTION_DEF ( char **alloca_matrix, (
int nrl,
int nrh,
int ncl,
int nch,
int entry_size,
char **m))
{
/* Allocates a matrix with the range [nrl .. nrh][ncl .. nch] */
int i;
char *start; /* Address of the [0][0] entry */
int rdim = nrh - nrl + 1; /* # entries in a row */
int cdim = nch - ncl + 1; /* # entries in a column */
/* Adjust the start of arrays, and adjust m */
start = (char *) (m + rdim); /* Count forward rdim (MTYPE *) */
start -= (ncl + nrl * cdim) * entry_size; /* Count back units of entry_size*/
m -= nrl;
/* Fill out the addresses */
for (i=nrl; i<=nrh; i++)
m[i] = &(start[i*cdim*entry_size]);
/* Return pointer to the pointer array */
return (m);
}
FUNCTION_DEF ( char **malloc_matrix, (
int nrl,
int nrh,
int ncl,
int nch,
int entry_size))
{
/* Allocates a matrix with the range [nrl .. nrh][ncl .. nch] */
int i;
char **m;
char *start; /* Address of the [0][0] entry */
int rdim = nrh - nrl + 1; /* # entries in a row */
int cdim = nch - ncl + 1; /* # entries in a column */
int nptrs = rdim;
int size;
/*
* We must make sure that double words are allocated on double word
* boundaries
*/
/* We must make sure that we have an even number of pointers */
if (entry_size > sizeof(char *) && nptrs % 2 != 0)
nptrs += 1;
/* Now, compute the size of space to be allocated */
size = nptrs*sizeof(char *) + rdim*cdim*entry_size;
/* Allocate m */
m = (char **) MALLOC(size);
start = (char *) (m + nptrs); /* Count forward nptrs (MTYPE *) */
start -= (ncl + nrl * cdim) * entry_size; /* Count back units of entry_size*/
m -= nrl;
/* Fill out the addresses */
for (i=nrl; i<=nrh; i++)
m[i] = &(start[i*cdim*entry_size]);
/* Return pointer to the pointer array */
return (m);
}
FUNCTION_DEF ( void fprint_intmatrix, (
char *format,
FILE *afile,
int **a,
int nrl,
int nrh,
int ncl,
int nch))
{
/* Print out a matrix of int values */
int i, j;
for (i=nrl; i<=nrh; i++)
{
for (j=ncl; j<=nch; j++)
fprintf (afile, format, a[i][j]);
fprintf (afile, "\n");
}
fprintf (afile, "\n");
}
FUNCTION_DEF ( void fprint_shortmatrix, (
char *format,
FILE *afile,
short **a,
int nrl,
int nrh,
int ncl,
int nch))
{
/* Print out a matrix of short values */
int i, j;
for (i=nrl; i<=nrh; i++)
{
for (j=ncl; j<=nch; j++)
fprintf (afile, format, a[i][j]);
fprintf (afile, "\n");
}
fprintf (afile, "\n");
}
FUNCTION_DEF ( void fprint_floatmatrix, (
char *format,
FILE *afile,
float **a,
int nrl,
int nrh,
int ncl,
int nch))
{
/* Print out a matrix of float values */
int i, j;
for (i=nrl; i<=nrh; i++)
{
for (j=ncl; j<=nch; j++)
fprintf (afile, format, a[i][j]);
fprintf (afile, "\n");
}
fprintf (afile, "\n");
}
FUNCTION_DEF ( void fprint_doublematrix, (
char *format,
FILE *afile,
double **a,
int nrl,
int nrh,
int ncl,
int nch))
{
/* Print out a matrix of double values */
int i, j;
for (i=nrl; i<=nrh; i++)
{
for (j=ncl; j<=nch; j++)
fprintf (afile, format, a[i][j]);
fprintf (afile, "\n");
}
fprintf (afile, "\n");
}
FUNCTION_DEF ( void fprint_matrix, (
char *format,
FILE *afile,
double **a,
int nrl,
int nrh,
int ncl,
int nch))
{
/* Print out a matrix of double values */
int i, j;
for (i=nrl; i<=nrh; i++)
{
for (j=ncl; j<=nch; j++)
fprintf (afile, format, a[i][j]);
fprintf (afile, "\n");
}
fprintf (afile, "\n");
}
FUNCTION_DEF ( void print_matrix, (
FILE *afile,
double **a,
int nrl,
int nrh,
int ncl,
int nch))
{
/* Print out a matrix of double values */
int i, j;
for (i=nrl; i<=nrh; i++)
{
int count = 0;
for (j=ncl; j<=nch; j++)
{
if (++count == 8)
{
fprintf (afile, "\n\t");
count = 1;
}
fprintf (afile,"%12.6f ", a[i][j]);
}
fprintf (afile, "\n");
}
fprintf (afile, "\n");
}
FUNCTION_DEF ( void fprint_vector, (
char *format,
FILE *afile,
double *v,
int nrl,
int nrh))
{
/* Print out a vector of double values */
int i;
for (i=nrl; i<=nrh; i++)
fprintf (afile, format, v[i]);
fprintf (afile, "\n");
}
FUNCTION_DEF ( void print_vector, (
FILE *afile,
double *v,
int nrl,
int nrh))
{
/* Print out a vector of double values */
int i;
int count = 0;
for (i=nrl; i<=nrh; i++)
{
if (++count == 8)
{
fprintf (afile, "\n\t");
count = 1;
}
fprintf (afile, "%9.6f ", v[i]);
}
fprintf (afile, "\n\n");
}
FUNCTION_DEF ( void matrix_prod, (
double **A,
double **B,
int l,
int m,
int n,
double **X))
{
/* Multiplies a k x m matrix (A) by an m x n matrix (B) to get matrix X */
/* All indices start at 1 */
/* The goal matrix may be the same as one of the operands */
int i, j, k;
double **T = MATRIX (1, l, 1, n, double); /* Temporary product */
/* Clear the matrix */
for (i=1; i<=l; i++)
for (j=1; j<=n; j++)
T[i][j] = 0.0;
/* Do the multiplication */
for (i=1; i<=l; i++)
for (j=1; j<=n; j++)
for (k=1; k<=m; k++)
T[i][j] += A[i][k] * B[k][j];
/* Clear the matrix */
for (i=1; i<=l; i++)
for (j=1; j<=n; j++)
X[i][j] = T[i][j];
/* Free T */
FREE (T, 1);
}
FUNCTION_DEF ( void vector_matrix_prod, (
double *b,
double **A,
int l,
int m,
double *x))
{
/* Multiplies a vector (b) by an l x m matrix (A) to get vector x */
/* All indices start at 1 */
/* The goal vector may be the same as b */
int i, j, k;
double *T = VECTOR (1, m, double); /* Temporary product */
/* Clear the vector */
for (i=1; i<=m; i++)
T[i] = 0.0;
/* Do the multiplication */
for (i=1; i<=m; i++)
for (k=1; k<=l; k++)
T[i] += b[k] * A[k][i];
/* Transfer the vector */
for (i=1; i<=m; i++)
x[i] = T[i];
/* Free T */
FREE (T, 1);
}
FUNCTION_DEF ( void matrix_vector_prod, (
double **A,
double *b,
int l,
int m,
double *x))
{
/* Multiplies a l x m matrix (A) by an m vector (b) to get vector x */
/* All indices start at 1 */
/* The goal vector may be the same as b */
int i, j, k;
double *T = VECTOR (1, l, double); /* Temporary product */
/* Clear the vector */
for (i=1; i<=l; i++)
T[i] = 0.0;
/* Do the multiplication */
for (i=1; i<=l; i++)
for (k=1; k<=m; k++)
T[i] += A[i][k] * b[k];
/* Transfer the vector */
for (i=1; i<=l; i++)
x[i] = T[i];
/* Free T */
FREE (T, 1);
}
FUNCTION_DEF ( double inner_product, (
double *A,
double *B,
int n))
{
/* Returns the inner product of two vectors */
int i;
double prod = 0.0;
/* Do the multiplication */
for (i=1; i<=n; i++)
prod += A[i] * B[i];
/* Return the product */
return (prod);
}
FUNCTION_DEF ( double vector_length, (
double *A,
int n))
{
/* Returns the length of the vector */
return (sqrt (inner_product (A, A, n)));
}
FUNCTION_DEF ( void transpose, (
double **A, /* The matrix to be transposed */
double **At, /* Its transpose (returned) */
int nrl, int nrh, /* X dimension */
int ncl, int nch /* Y dimension */
))
{
/* Does the transpose of a matrix. */
/* Goal and origin may be the same. */
int i, j;
double **temp;
int del = 0;
/* First does the diagonal entries */
if (A == At)
{
temp = MATRIX (nrl, nrh, ncl, nch, double);
del = 1;
for (i=nrl; i<=nrh; i++) for (j=1; j<=nch; j++)
temp[i][j] = A[i][j];
}
else temp = A;
/* Now swaps the off-diagonal entries */
for (i=nrl; i<=nrh; i++) for (j=ncl; j<=nch; j++)
At[j][i] = temp[i][j];
/* Free the matrix if required */
if(del) FREE(temp, nrl);
}
FUNCTION_DEF ( int matrix_inverse, (
double **A, /* The matrix to be inverted */
double **Ainv, /* The inverse computed */
int n /* Dimension of the matrix */
))
{
/* Inverts the matrix */
/* Returns 0 if matrix is not invertible */
int i, j;
/* First copy the matrix to Ainv */
for (i=1; i<=n; i++) for (j=1; j<=n; j++)
Ainv[i][j] = A[i][j];
/* Now invert it on the spot */
return gaussj2_check (Ainv, n, (double **)0, 0);
}
FUNCTION_DEF ( char **sub_matrix, (
char **A, /* Base array */
int nrl, int ncl, /* Origin of the new array */
int rdim, /* Number of rows of the sub-array */
int orgi, int orgj, /* Starting index for new array */
int entry_size /* Size of an entry */
))
{
/* Allocates a sub-matrix of A with the given range */
int i;
char **m;
int offset;
int size = rdim*sizeof(char *);
/* Allocate m */
m = (char **) MALLOC(size);
m -= nrl ;
/* Fill out the addresses */
offset = (ncl-orgj) * (entry_size/sizeof(char));
for (i=nrl; i<nrl+rdim; i++)
m[i] = A[i] + offset;
/* Return pointer to the pointer array */
return (m+nrl-orgi);
}
FUNCTION_DEF ( double det3x3, (
double **A))
{
/* Takes a 3 by 3 determinant */
return A[1][1]*A[2][2]*A[3][3] + A[1][2]*A[2][3]*A[3][1] +
A[1][3]*A[2][1]*A[3][2] - A[1][1]*A[2][3]*A[3][2] -
A[1][2]*A[2][1]*A[3][3] - A[1][3]*A[2][2]*A[3][1];
}
FUNCTION_DEF ( void cofactor_matrix_3x3, (
double **A,
double **Aadj))
{
/* Takes the cofactor matrix of a given 3x3 matrix */
int i, j;
/* First, take a copy of A */
double **A2 = MATRIX (1, 3, 1, 3, double);
for (i=1; i<=3; i++) for (j=1; j<=3; j++)
A2[i][j] = A[i][j];
/* Now take the cofactors */
Aadj[1][1] = A2[2][2]*A2[3][3] - A2[2][3]*A2[3][2];
Aadj[1][2] = A2[2][3]*A2[3][1] - A2[2][1]*A2[3][3];
Aadj[1][3] = A2[2][1]*A2[3][2] - A2[2][2]*A2[3][1];
Aadj[2][1] = A2[3][2]*A2[1][3] - A2[3][3]*A2[1][2];
Aadj[2][2] = A2[3][3]*A2[1][1] - A2[3][1]*A2[1][3];
Aadj[2][3] = A2[3][1]*A2[1][2] - A2[3][2]*A2[1][1];
Aadj[3][1] = A2[1][2]*A2[2][3] - A2[1][3]*A2[2][2];
Aadj[3][2] = A2[1][3]*A2[2][1] - A2[1][1]*A2[2][3];
Aadj[3][3] = A2[1][1]*A2[2][2] - A2[1][2]*A2[2][1];
/* Free the temporary matrix */
FREE (A2, 1);
}
FUNCTION_DEF ( void matrix_copy, (
double **A, double **B, /* The matrices to copy from and to */
int i0, int i1, /* Row bounds */
int j0, int j1 /* Column bounds */
))
{
/* Copies a matrix */
int i, j;
for (i=i0; i<=i1; i++)
for (j=j0; j<=j1; j++)
B[i][j] = A[i][j];
}
FUNCTION_DEF ( void cross_product, (
double *A, double *B, /* Two vectors 1..3 */
double *C /* Their cross product */
))
{
C[1] = A[2]*B[3] - A[3]*B[2];
C[2] = A[3]*B[1] - A[1]*B[3];
C[3] = A[1]*B[2] - A[2]*B[1];
}
FUNCTION_DEF ( int lin_solve, (
double **A, /* The matrix of coefficients */
double *x, /* The vector of unknowns */
double *b, /* The goal vector */
int n /* Size of the system */
))
{
/* Solves a square linear system of equations */
int i, j;
int return_val;
double **A2 = MATRIX (1, n, 1, n, double);
double *b2 = VECTOR (1, n, double);
/* Copy the values */
for (i=1; i<=n; i++) for (j=1; j<=n; j++)
A2[i][j] = A[i][j];
for (i=1; i<=n; i++)
b2[i] = b[i];
/* Now solve using gaussj */
return_val = gaussj_check (A2, n, b2);
/* Now copy the values */
if (return_val)
for (i=1; i<=n; i++)
x[i] = b2[i];
/* Free the temporary matrices */
FREE (A2, 1);
FREE (b2, 1);
/* Return the appropriate value */
return return_val;
}
FUNCTION_DEF ( void identity_matrix, (
double **A,
int dim))
{
/* Initializes the matrix A to be the identity */
int i, j;
for (i=1; i<=dim; i++)
{
A[i][i] = 1.0;
for (j=i+1; j<=dim; j++)
A[i][j] = A[j][i] = 0.0;
}
}