mops.h Source File

00001 //
00002 // mops.h --- block matrix operations
00003 //
00004 // Copyright (C) 1997 Limit Point Systems, Inc.
00005 //
00006 // Author: Edward Seidl <seidl@janed.com>
00007 // Maintainer: LPS
00008 //
00009 // This file is part of the SC Toolkit.
00010 //
00011 // The SC Toolkit is free software; you can redistribute it and/or modify
00012 // it under the terms of the GNU Library General Public License as published by
00013 // the Free Software Foundation; either version 2, or (at your option)
00014 // any later version.
00015 //
00016 // The SC Toolkit is distributed in the hope that it will be useful,
00017 // but WITHOUT ANY WARRANTY; without even the implied warranty of
00018 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00019 // GNU Library General Public License for more details.
00020 //
00021 // You should have received a copy of the GNU Library General Public License
00022 // along with the SC Toolkit; see the file COPYING.LIB.  If not, write to
00023 // the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
00024 //
00025 // The U.S. Government is granted a limited license as per AL 91-7.
00026 //
00027 
00028 #ifndef _math_scmat_mops_h
00029 #define _math_scmat_mops_h
00030 
00031 #define D1 32
00032 
00033 // copy a chunk of rectangular matrix source into dest.  dest is D1xD1, and is
00034 // padded with zeros
00035 
00036 static inline void
00037 copy_block(double **dest, double **source,
00038            int istart, int ni, int jstart, int nj)
00039 {
00040   int ii,jj;
00041   
00042   for (ii=0; ii < ni; ii++) {
00043     double *di = dest[ii];
00044     double *si = &source[istart+ii][jstart];
00045     for (jj=0; jj < nj; jj++)
00046       di[jj] = si[jj];
00047     for (; jj < D1; jj++)
00048       di[jj] = 0;
00049   }
00050 
00051   int left=D1-ii;
00052   if (left)
00053     memset(dest[ii], 0, sizeof(double)*left*D1);
00054 }
00055 
00056 static inline void
00057 copy_trans_block(double **dest, double **source,
00058                  int istart, int ni, int jstart, int nj)
00059 {
00060   int ii,jj;
00061   
00062   memset(dest[0], 0, sizeof(double)*D1*D1);
00063 
00064   for (jj=0; jj < nj; jj++) {
00065     double *sj = &source[jstart+jj][istart];
00066     for (ii=0; ii < ni; ii++)
00067       dest[ii][jj] = sj[ii];
00068   }
00069 }
00070 
00071 // copy a chunk of symmetric matrix source into dest.  dest is D1xD1, and is
00072 // padded with zeros
00073 static inline void
00074 copy_sym_block(double **dest, double **source,
00075                int istart, int ni, int jstart, int nj)
00076 {
00077   int ii,jj;
00078 
00079   for (ii=0; ii < ni; ii++) {
00080     double *di = dest[ii];
00081     double *si = &source[istart+ii][jstart];
00082     
00083     if (jstart < istart)
00084       for (jj=0; jj < nj; jj++)
00085         di[jj] = si[jj];
00086     else if (jstart==istart)
00087       for (jj=0; jj <= ii; jj++)
00088         di[jj] = dest[jj][ii] = si[jj];
00089     else
00090       for (jj=0; jj < nj; jj++)
00091         di[jj] = source[jstart+jj][istart+ii];
00092 
00093     for (jj=nj; jj < D1; jj++)
00094       di[jj] = 0;
00095   }
00096 
00097   int left=D1-ii;
00098   if (left)
00099     memset(dest[ii], 0, sizeof(double)*left*D1);
00100 }
00101 
00102 static inline void
00103 return_block(double **dest, double **source,
00104              int istart, int ni, int jstart, int nj)
00105 {
00106   int ii,jj;
00107 
00108   for (ii=0; ii < ni; ii++)
00109     for (jj=0; jj < nj; jj++)
00110       dest[istart+ii][jstart+jj] = source[ii][jj];
00111 }
00112 
00113 // a, b, and c are all D1xD1 blocks
00114 static inline void
00115 mult_block(double **a, double **b, double **c, int ni, int nj, int nk)
00116 {
00117   int ii,jj,kk;
00118   double t00,t10,t20,t30;
00119   double *a0, *a1, *a2, *a3;
00120   double *c0, *c1, *c2, *c3;
00121 
00122   for (ii=0; ii < ni; ii += 4) {
00123     a0=a[ii]; a1=a[ii+1]; a2=a[ii+2]; a3=a[ii+3];
00124     c0=c[ii]; c1=c[ii+1]; c2=c[ii+2]; c3=c[ii+3];
00125 
00126     for (jj=0; jj < nj; jj++) {
00127       double *bt = b[jj];
00128       t00=c0[jj]; t10=c1[jj]; t20=c2[jj]; t30=c3[jj];
00129 
00130       for (kk=0; kk < nk; kk += 2) {
00131         register double b0=bt[kk], b1=bt[kk+1];
00132         t00 += a0[kk]*b0 + a0[kk+1]*b1;
00133         t10 += a1[kk]*b0 + a1[kk+1]*b1;
00134         t20 += a2[kk]*b0 + a2[kk+1]*b1;
00135         t30 += a3[kk]*b0 + a3[kk+1]*b1;
00136       }
00137 
00138       c0[jj]=t00;
00139       c1[jj]=t10;
00140       c2[jj]=t20;
00141       c3[jj]=t30;
00142     }
00143   }
00144 }
00145 
00146 #endif
00147 
00148 // Local Variables:
00149 // mode: c++
00150 // c-file-style: "ETS"
00151 // End: