---

CTBdot.cxx

Go to the documentation of this file.

/*----------------------------------------------------------------------------*/
/* C Tool Box: Designed and implemented by:                                   */
/*    Walter F.J. Mueller   Gesellschaft fuer Schwerionenforschung (GSI)      */
/*                          Planckstrasse 1, D-64291 Darmstadt, Germany       */
/*                  Email:  W.F.J.Mueller@gsi.de                              */
/*                  WWW:    http://www-kp3.gsi.de/www/kp3/people/mueller.html */
/*------------------------------------------+---------------------------------*/

#include "CTB.hxx"
#include "CTBtraits.hxx"
#include "CTBnum.hxx"
#include "CTBexceptionConformity.hxx"
#include "CTBdot.hxx"

//------------------------------------------+-----------------------------------
template <class T>
void CTBdot(CTBvector<T>& r, const CTBmatrix<T>& m, const CTBvector<T>& v)
{
  CTBint   i_nr = m.NRow();                 // #rows of m
  CTBint   i_nc = m.NColumn();              // #cols of m
  const T* p_m  = m;                        // running ptr in m
  T*       p_r;                             // ptr to result vector

  if (i_nc != v.Size()) 
    throw CTBexceptionConformity(i_nc, v.Size(),
                                 "CTBdot(CTBvector&,CTBmatrix&,CTBvector&");
  
  if (r.Size() != i_nr) {                   // resize output if needed
    r.Resize(0);                            // zero first to avoid copy
    r.Resize(i_nr);
  }
  
  p_r  = r;

  for (CTBint i = 0; i < i_nr; i++) {       // loop over rows
    const T* p_v = v;
    CTBtraits<T>::Sum_t f_sum = 0.;

    for (CTBint j = 0; j < i_nc; j++) f_sum += *p_m++ * *p_v++;
    *p_r++ = f_sum;
  }
  return;
}

//------------------------------------------+-----------------------------------
template <class T>
void CTBdot(CTBmatrix<T>& r, const CTBmatrix<T>& m1, const CTBmatrix<T>& m2)
{
  CTBint   i_nr = m1.NRow();                // #row of result
  CTBint   i_nc = m2.NColumn();             // #col of result
  CTBint   i_nk = m1.NColumn();             // length of dot products
  const T* p_m2 = m2;                       // ptr to current column in m2
  T*       p_r;                             // ptr to current column in r

  if (i_nk != m2.NRow()) 
    throw CTBexceptionConformity(i_nk, m2.NRow(),
                                 "CTBdot(CTBmatrix&,CTBmatrix&,CTBmatrix&");
  
  if (r.NRow() != i_nr || r.NColumn() != i_nc) {
    r.Resize(0,0);
    r.Resize(i_nr,i_nc);
  }

  p_r  = r;

  for (CTBint i = 0; i < i_nr; i++) {
    const T* p_m1 = m1;                     // running ptr in m1
    T*       p_vr = p_r;                    // running ptr to column in r

    for (CTBint j = 0; j < i_nc; j++) {
      const T* p_v  = p_m2;                 // running ptr to column in m2
      CTBtraits<T>::Sum_t f_sum = 0.;

      for (CTBint k = 0; k < i_nk; k++) {
        f_sum += *p_m1++ * *p_v;
        p_v   += i_nc;                      // to element in m2 column
      }

      *p_vr = f_sum;
      p_vr += i_nc;
    }

    p_m2++;                                 // point to next column in m2
    p_r++;                                  // point to next column in r
  }
  return;
}

---