imat10.cpp

Go to the documentation of this file.

#include "fvar.hpp"

ivector diagonal(const imatrix& matrix)
{
  int min = matrix.rowmin();
  int max = matrix.rowmax();
  if (min != matrix.colmin() || max != matrix.colmax())
  {
    cerr << "Error: imatrix is not square"
         << " in diagonal(const imatrix&).\n";
    ad_exit(1);
  }
  ivector vector(min, max);
  for (int i = min; i <= max; ++i)
  {
    vector(i) = matrix(i,i);
  }
  return vector;
}
int operator*(const ivector& t1, const ivector& t2)
{
  if (t1.indexmin() != t2.indexmin() ||  t1.indexmax() != t2.indexmax())
  {
    cerr << "Index bounds do not match"
         << " in operator*(const ivector&, const ivector&)\n";
    ad_exit(1);
  }
  int tmp = 0;
#ifdef OPT_LIB
  const int* pt1=&t1[t1.indexmin()];
  const int* pt1m=&t1[t1.indexmax()];
  const int* pt2=&t2[t2.indexmin()];
  do
  {
    tmp+=*pt1++ * *pt2++;
  }
  while (pt1<=pt1m);

#else
  #ifndef USE_ASSEMBLER
    int min=t1.indexmin();
    int max=t1.indexmax();
    for (int i=min; i<=max; i++)
    {
      tmp+=t1[i]*t2[i];
    }
  #else
    int min=t1.indexmin();
    int n=t1.indexmax()-min+1;
    dp_dotproduct(&tmp,&(t1(min)),&(t2(min)),n);
  #endif
#endif

  return tmp;
}
imatrix operator*(const imatrix& a, const imatrix& b)
{
  int min = a.rowmin();
  int max = a.rowmax();
  if (min != b.colmin() || max != b.colmax())
  {
    cerr << " Incompatible array bounds"
         << " in operator*(const imatrix&, const imatrix&)\n";
    ad_exit(1);
  }
  imatrix results(min, max, min, max);
  for (int j = min; j <= max; ++j)
  {
    ivector col = column(b, j);
    for (int i = min; i <= max; ++i)
    {
      results(i, j) = a(i) * col;
    }
  }
  return results;
}