/* Copyright (C) 1982-1999 Ken Turkowski.
*
* All rights reserved.
*
* Warranty Information
* Even though I have reviewed this software, I makes no warranty
* or representation, either express or implied, with respect to this
* software, its quality, accuracy, merchantability, or fitness for a
* particular purpose. As a result, this software is provided "as is,"
* and you, its user, are assuming the entire risk as to its quality
* and accuracy.
*
* This code may be used and freely distributed as long as it includes
* this copyright notice and the above warranty information.
*
********************************************************************************
* LinearTransform
* General-purpose linear transformation utility
*
* Examples:
* LinearTransform(v, v, d, 1, 3, 1); // Squared norm of a 3-vector, yielding a scalar.
* LinearTransform(u, v, d, 1, 3, 1); // Dot product of two 3-vectors, yielding a scalar.
* LinearTransform(u, v, d, 3, 1, 3); // Tensor product of two 3-vectors, yielding a 3x3 2-rank tensor.
* LinearTransform(v, M, w, 1, 3, 3); // Transforms a 1x3 row vector by a 3x3 matrix, yielding a 1x3 row vector.
* LinearTransform(V, M, W, 16, 4, 4); // Transforms 16 1x4 row vectors by a 4x4 matrix, yielding 16 1x4 row vectors.
* LinearTransform(V, T, U, 3, 4, 3); // Multiplies a 3x4 matrix by a 4x3 matrix yielding a 3x3 matrix.
********************************************************************************/
#include
#define real double /* The data types of L, R and P */
/* Prototype */
void LinearTransform(real *L, real *R, real *P, unsigned int nRow, unsigned int lCol, unsigned int rCol);
/********************************************************************************
* LinearTransform multiplies an (nRow x lCol) matrix by a (lCol x rCol) matrix,
* yielding an (nRow x rCol) matrix.
* The source matrices are given in L and R, and the result is returned in P.
********************************************************************************/
void
LinearTransform(
real *L, /* nRow x lCol */
real *R, /* lCol x rCol */
register real *P, /* nRow x rCol */
register unsigned int nRow, /* Number of rows in L and P */
unsigned int lCol, /* Number of columns in L and number of rows in R */
unsigned int rCol /* Number of columns in R and P */
)
{
register real *ap, *bp;
register unsigned int k, j, i;
long double sum; /* Extended precision for intermediate results */
for (i = nRow; i--; L += lCol) { /* Each row in L */
for (j = 0; j < rCol; j++) { /* Each column in R */
ap = L; /* Left of ith row of L */
bp = R + j; /* Top of jth column of R */
sum = 0.0;
for (k = lCol; k--; bp += rCol)
sum += *ap++ * (*bp); /* *P += L[i'][k'] * R[k'][j]; */
*P++ = sum;
}
}
}