imat10.cpp

Go to the documentation of this file.

/*
 * $Id$
 *
 * Author: David Fournier
 * Copyright (c) 2008-2012 Regents of the University of California
 */
#include "fvar.hpp"

ivector diagonal(const imatrix& m)
{
  if (m.indexmin() != m.colmin() || m.indexmax() != m.colmax())
  {
    cerr << "Error matrix not square in function diagonal" << endl;
    exit(21);
  }
  int mmin=m.indexmin();
  int mmax=m.indexmax();
  ivector tmp(mmin,mmax);
  for (int i=mmin;i<=mmax;i++)
    tmp(i)=m(i,i);
  return tmp;
}
int operator*(const ivector& t1, const ivector& t2)
{
  if (t1.indexmin() != t2.indexmin() ||  t1.indexmax() != t2.indexmax())
  {
    cerr << "Index bounds do not match in "
      "ivector 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& m1, const imatrix& m2 )
{
  if (m1.colmin() != m2.rowmin() || m1.colmax() != m2.rowmax())
  {
    cerr << " Incompatible array bounds in "
      "imatrix  operator * (const imatrix& x, const imatrix& m)\n";
    ad_exit(21);
  }
  imatrix tmp(m1.rowmin(),m1.rowmax(), m2.colmin(), m2.colmax());
  for (int j=m2.colmin(); j<=m2.colmax(); j++)
  {
    ivector m2col=column(m2,j);
    for (int i=m1.rowmin(); i<=m1.rowmax(); i++)
    {
      tmp(i,j) = m1(i) * m2col;
    }
  }
  return tmp;
}