/* mat.c Simple matrix manipulation routines.
* © Copyright 1999 by John Halleck
* All rights reserved.
*/
/* Version of August 18th, 1999 */
/* Initial Test Version
* 0.11 - John Halleck, September 6th
* added matrix subtract.
* 0.10 - John Halleck, July 26th, 1999
* removed embed and extract, and moved them to a partitioning package.
* 0.09 - John Halleck, July 24th, 1999
* Added embed and extract operations to support matrix partitioning.
* 0.08 - John Halleck, July 17, 1999
* Stripped out most of the guts. Pulled vector stuff.
* Only the stuff we use is left. That should make
* it less work to make matching symmetric versions.
* 0.07 - John Halleck, July 16, 1999
* Stripped some functions out.
* 0.06 - John Halleck, June 30th, 1999
* Regularized form to match rest of package.
* 0.05 - John Halleck, June 29th, 1999
* pulled various max functions. Added copy transpose.
* 0.04 - John Halleck, June 21st, 1999
* added rowmuls
* 0.03 - John Halleck June 8th, 1999
* added some routines for interface completeness.
* 0.02 - John Halleck June 5th, 1999
* added row vector vs col vector stuff.
* 0.01 - John Halleck June 4, 1999
* Initial creation.
*/
#include "errors.h"
/* Standard error codes. */
#include "mat.h"
/* Matrix definitions. */
/* ========================= Set to a constant ============================== */
/* ------------- matzero ---------------- */
/* Make a zero matrix. */
error matzero (int rows, int cols, matrix result) {
int size; double *out;
if (!result) return ERRnil; /* We have to have a matrix */
if (rows<1 || cols<1) return ERRsize; /* Of reasonable size */
size = rows*cols; out = (double *) result;
while (size--) *out++ = 0.0;
return NoError;
}
/* ------------- matidn ----------------- */
/* matrix identity */
error matidn (int size, matrix result)
{
double *lmatrix;
int index; long fullsize;
if (!result) return ERRnil; /* Must have an array */
if (size < 1) return ERRsize; /* With an actual size */
fullsize = size*size;
index = fullsize; lmatrix = (double *) result;
while (index--) *lmatrix++ = 0.0;
index = size; lmatrix = (double *) result;
while (index--) { *lmatrix = 1.0; lmatrix += size + 1; }
return NoError;
}
/* ========================= Copy =========================================== */
/* ------------- matcpy ----------------- */
/* Copy one matrix to another. */
error matcpy (int rows, int cols, matrix result, matrix source) {
long size; double *final, *initial;
if (!result || !source) return ERRnil; /* Must have arrays */
if (rows < 1 || cols < 1) return ERRsize; /* of reasonable size */
size = rows*cols; final = (double *)result; initial = (double *) source;
while (size--) { *final++ = *initial++; }
return NoError;
}
/* ---------- Copy Transpose ------------ */
/* Copy a matrix to its transpose */
/* Source and destination must be distinct. */
/* rows and columns refer to destination. */
error matcpyt (int rows, int cols, matrix result, matrix source) {
double *final, *initial;
int i, j;
if (!result || !source) return ERRnil; /* Must have arrays */
if (rows < 1 || cols < 1) return ERRsize; /* of reasonable size */
if (source == result) return ERRsame; /* And they must be distinct */
final = (double *)result; initial = (double *) source;
for (i=0;i<rows;i++) for (j=0;j<cols;j++) {
*final++ = *(initial + j*rows + i);
}
return NoError;
}
/* ======================== Multiply by scalar ============================== */
/* ------------------ matmuls ---------------- */
/* Multiply matrix by scalar. */
error matscl (int rows, int cols, matrix result, matrix given, double scalar) {
long size; double *dst, *src;
if (!result || !given) return ERRnil; /* Must have an array */
if (rows < 1 || cols < 1) return ERRsize; /* of reasonable size */
dst = (double *) result; src = (double *) given;
size = rows * cols;
while (size--) *dst++ = scalar * *src++;
if (scalar == 0.0) return ERRmeaning;
return NoError;
}
/* ========================== Add =========================================== */
/* result = A + B */
/* ------------------ matadd ------------------ */
/* Add two matrices. */
error matadd (int rows, int cols, matrix result, matrix A, matrix B) {
long size; double *final, *locala, *localb;
if (!result || !A || !B) return ERRnil; /* Must have arrays */
if (rows < 1 || cols < 1) return ERRsize; /* of reasonable size */
final = (double *) result; locala = (double *) A; localb = (double *) B;
size = rows * cols;
while (size--) *final++ = *locala++ + *localb++;
return NoError;
}
/* result = A - B */
/* ----------------- matsub ------------------ */
error matsub (int rows, int cols, matrix result, matrix A, matrix B) {
long size; double *final, *locala, *localb;
if (!result || !A || !B) return ERRnil; /* Must have arrays */
if (rows < 1 || cols < 1) return ERRsize; /* of reasonable size */
final = (double *) result; locala = (double *) A; localb = (double *) B;
size = rows * cols;
while (size--) *final++ = *locala++ - *localb++;
return NoError;
}
/* ========================== Multiply ===================================== */
/* ------------------ matmul ----------------- */
/* Multiply matrices. result = matrix a multiplied by matrix b; */
/* result MUST be a matrix distinct from either A or B. */
error matmul (int rrows, int rcols, matrix result,
int isize, matrix A, matrix B)
/* result is rrows x rcols */
/* A is rrows x isize */
/* B is isize x rcols */
{
double mterm; /* Term we are collecting */
double *rowpoint, *colpoint; /* pointer to start of current row and column */
double *ipointa, *ipointb; /* current locations for inner product. */
double *final;
int thisrow, thiscol, index;
if (!result || !A || !B) return ERRnil; /* Must have arrays */
if (rrows < 1 || rcols < 1 || isize < 1)
return ERRsize; /* With real size */
if (result == A || result == B) return ERRsame; /* and result is distinct */
final = (double *) result; rowpoint = (double *) A;
for (thisrow = 0; thisrow < rrows; thisrow++) {
colpoint = (double *) B;
for (thiscol = 0; thiscol < rcols; thiscol++) {
mterm = 0.0;
ipointa = rowpoint; ipointb = colpoint;
for (index=0; index < isize; index++) {
mterm += *ipointa++ * *ipointb; ipointb += rcols;
}
*final++ = mterm;
colpoint++; /* On to next column */
}
rowpoint += isize; /* Pointer to first item of this row in A */
}
return NoError;
}
/* ------------------ matmult ----------------- */
/* Multiply matrix by a transpose. */
error matmult (int rrows, int rcols, matrix result,
int isize, matrix A, matrix B)
/* result is rrows x rcols */
/* A is rrows x isize */
/* B is rcols x isize <= */
{
double mterm; /* Term we are collecting */
double *rowpoint, *colpoint; /* pointer to start of current row and column */
double *ipointa, *ipointb; /* current locations for inner product. */
double *final;
int thisrow, thiscol, index;
if (!result || !A || !B) return ERRnil; /* Must have arrays */
if (rrows < 1 || rcols < 1 || isize < 1)
return ERRsize; /* With real size */
if (result == A || result == B) return ERRsame; /* and result is distinct */
final = (double *) result;
rowpoint = (double *) A; colpoint = (double *) B;
for (thisrow = 0; thisrow < rrows; thisrow++) {
for (thiscol = 0; thiscol < rcols; thiscol++) {
mterm = 0.0;
ipointa = rowpoint; ipointb = colpoint;
for (index=0; index < isize; index++) {
mterm += *ipointa++ * *ipointb++;
}
*final++ = mterm;
colpoint += isize;
}
rowpoint += isize; /* Pointer to first item of this row in A */
colpoint = (double *) B;
}
return NoError;
}
This page is http://www.cc.utah.edu/~nahaj/cave/survey/code/c/mat.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