Main Page | Namespace List | Class Hierarchy | Alphabetical List | Class List | Directories | File List | Namespace Members | Class Members | File Members | Related Pages

matrix2.cpp

Go to the documentation of this file.
00001 // This may look like C code, but it is really -*- C++ -*-
00002 /*
00003  ************************************************************************
00004  *
00005  *                      Linear Algebra Package
00006  *
00007  *              Basic linear algebra operations, level 2
00008  *             Matrix transformations and multiplications
00009  *                         of various types
00010  *
00011  * $Id: matrix2.cpp,v 1.1 2004/05/21 21:02:52 maxx Exp $
00012  *
00013  ************************************************************************
00014  */
00015 
00016 #include "LinAlg.h"
00017 #include <math.h>
00018 #if defined(WIN32)
00019 #include <malloc.h>
00020 #else
00021 #include <alloca.h>
00022 #endif
00023 
00024 namespace linalg 
00025 {
00026     using namespace linalg;
00027     
00028 /*
00029  *------------------------------------------------------------------------
00030  *                              Transpositions
00031  */
00032                                 // Transpose a matrix
00033 void LazyTransposedMatrix::fill_in(Matrix& m) const
00034 {
00035   register const REAL * rsp = proto.elements;   // Row source pointer
00036   register REAL * tp = m.elements;
00037 
00038                                 // Matrix m is traversed in the
00039                                 // natural, column-wise way, whilst the source
00040                                 // (prototype) matrix is scanned row-by-row
00041   while( tp < m.elements + m.nelems )
00042   {
00043     register const REAL * sp = rsp++;   // sp = @proto[j,i] for i=0
00044                                         // Move tp to the next elem in the col,
00045     while( sp < proto.elements + proto.nelems )
00046        *tp++ = *sp, sp += proto.nrows; // sp to the next el in the curr row
00047   }
00048   assert( tp == m.elements + m.nelems && 
00049           rsp == proto.elements + proto.nrows );
00050 
00051 }
00052 
00053 /*
00054  *------------------------------------------------------------------------
00055  *                      General matrix multiplications
00056  */
00057 
00058                         // Compute target = target * source inplace.
00059                         // Strictly speaking, it can't be done inplace,
00060                         // though only one row of the target matrix needs
00061                         // to be saved.
00062                         // "Inplace" multiplication is only possible
00063                         // when the 'source' matrix is square
00064 Matrix& Matrix::operator *= (const Matrix& source)
00065 {
00066   is_valid();
00067   source.is_valid();
00068 
00069   if( row_lwb != source.col_lwb || ncols != source.nrows ||
00070       col_lwb != source.col_lwb || ncols != source.ncols )
00071     info(), source.info(),
00072     _error("matrices above are unsuitable for the inplace multiplication");
00073 
00074                                         // One row of the old_target matrix
00075 #ifdef WIN32
00076   REAL * const one_row = (REAL *)malloc(ncols*sizeof(REAL));
00077 #else
00078   REAL * const one_row = (REAL *)alloca(ncols*sizeof(REAL));
00079 #endif
00080   const REAL * one_row_end = &one_row[ncols];
00081 
00082   register REAL * trp = elements;       // Pointer to the i-th row
00083   for(; trp < &elements[nrows]; trp++)  // Go row-by-row in the target
00084   {
00085     register REAL *wrp, *orp;                   // work row pointers
00086     for(wrp=trp,orp=one_row; orp < one_row_end;)
00087       *orp++ = *wrp, wrp += nrows;              // Copy a row of old_target
00088 
00089     register REAL *scp=source.elements;         // source column pointer
00090     for(wrp=trp; wrp < elements+nelems; wrp += nrows)
00091     {
00092       register double sum = 0;                  // Multiply a row of old_target
00093       for(orp=one_row; orp < one_row_end;)      // by each col of source
00094         sum += *orp++ * *scp++;                 // to get a row of new_target
00095       *wrp = sum;
00096     }
00097   }
00098 
00099   return *this;
00100 }
00101 
00102                         // Compute C = A*B
00103                         // The matrix C must be already
00104                         // allocated, and it is *this
00105 void Matrix::mult(const Matrix& A, const Matrix& B)
00106 {
00107   A.is_valid();
00108   B.is_valid();
00109   is_valid();
00110   
00111   if( A.ncols != B.nrows || A.col_lwb != B.row_lwb )
00112     A.info(), B.info(),
00113     _error("matrices above cannot be multiplied");
00114   if( nrows != A.nrows || ncols != B.ncols ||
00115       row_lwb != A.row_lwb || col_lwb != B.col_lwb )
00116     A.info(),B.info(),info(),
00117     _error("product A*B is incompatible with the given matrix");
00118 
00119   register REAL * arp;                  // Pointer to the i-th row of A
00120            REAL * bcp = B.elements;     // Pointer to the j-th col of B
00121   register REAL * cp = elements;        // C is to be traversed in the natural
00122   while( cp < elements + nelems )       // order, col-after-col
00123   {
00124     for(arp = A.elements; arp < A.elements + A.nrows; )
00125     {
00126       register double cij = 0;
00127       register REAL * bccp = bcp;               // To scan the jth col of B
00128       while( arp < A.elements + A.nelems )      // Scan the i-th row of A and
00129         cij += *bccp++ * *arp, arp += A.nrows;  // the j-th col of B
00130       *cp++ = cij;
00131       arp -= A.nelems - 1;                      // arp points to (i+1)-th row
00132     }
00133     bcp += B.nrows;                     // We're done with j-th col of both
00134   }                                     // B and C. Set bcp to the (j+1)-th col
00135 
00136   assert( cp == elements + nelems && bcp == B.elements + B.nelems );
00137 }
00138 
00139 #if 0
00140                         // Create a matrix C such that C = A' * B
00141                         // In other words,
00142                         // c[i,j] = SUM{ a[k,i] * b[k,j] }
00143                         // Note, matrix C needs to be allocated
00144 void Matrix::_AtmultB(const Matrix& A, const Matrix& B)
00145 {
00146   A.is_valid();
00147   B.is_valid();
00148 
00149   if( A.nrows != B.nrows || A.row_lwb != B.row_lwb )
00150     A.info(), B.info(),
00151     _error("matrices above are unsuitable for A'B multiplication");
00152 
00153   allocate(A.ncols,B.ncols,A.col_lwb,B.col_lwb);
00154 
00155   register REAL * acp;                  // Pointer to the i-th col of A
00156            REAL * bcp = B.elements;     // Pointer to the j-th col of B
00157   register REAL * cp = elements;        // C is to be traversed in the natural
00158   while( cp < elements + nelems )       // order, col-after-col
00159   {
00160     for(acp = A.elements; acp<A.elements + A.nelems;)
00161     {                                   // Scan all cols of A
00162       register double cij = 0;                  
00163       register REAL * bccp = bcp;               // To scan the jth col of B
00164       for(register int i=0; i<A.nrows; i++)     // Scan the i-th row of A and
00165         cij += *bccp++ * *acp++;                        // the j-th col of B
00166       *cp++ = cij;
00167     }
00168     bcp += B.nrows;                     // We're done with j-th col of both
00169   }                                     // B and C. Set bcp to the (j+1)-th col
00170 
00171   assert( cp == elements + nelems && bcp == B.elements + B.nelems );
00172 }
00173 #endif
00174 }

Generated on Wed Dec 15 21:20:29 2004 for vuVolume by  doxygen 1.3.9.1