ADMB Documentation  11.1x.2735
 All Classes Files Functions Variables Typedefs Friends Defines
dmat1.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id: dmat1.cpp 2287 2014-09-04 18:12:38Z johnoel $
00003  *
00004  * Author: David Fournier
00005  * Copyright (c) 2008-2012 Regents of the University of California
00006  */
00011 #include "fvar.hpp"
00012 
00017 dvector operator*(const dvector& x, const dmatrix& m)
00018  {
00019 #ifdef DIAG
00020      if( heapcheck() == _HEAPCORRUPT )
00021      {
00022         if (ad_printf)
00023          (*ad_printf)( "Entering dvector * dvec dmat is corrupted.\n" );
00024      }
00025      else
00026      {
00027         if (ad_printf)
00028           (*ad_printf)( "Entering dvector * dvec dmat  Heap is OK.\n" );
00029      }
00030 #endif
00031 
00032    if (x.indexmin() != m.rowmin() || x.indexmax() != m.rowmax())
00033    {
00034      cerr << " Incompatible array bounds in "
00035      "dvector  operator * (const dvector& x,const dmatrix& m)\n";
00036      ad_exit(21);
00037    }
00038    dvector tmp(m.colmin(),m.colmax());
00039 
00040    for (int j=m.colmin(); j<=m.colmax(); j++)
00041    {
00042      tmp[j]=0;
00043      for (int i=x.indexmin(); i<=x.indexmax(); i++)
00044      {
00045        tmp[j]+=x[i]*m[i][j];
00046      }
00047    }
00048 #ifdef DIAG
00049      if( heapcheck() == _HEAPCORRUPT )
00050      {
00051         if (ad_printf)
00052           (*ad_printf)( "Leaving dvector * dvec dmat is corrupted.\n" );
00053      }
00054      else
00055      {
00056         if (ad_printf)
00057           (*ad_printf)( "Leaving dvector * dvec dmat  Heap is OK.\n" );
00058      }
00059 #endif
00060    return(tmp);
00061  }
00062 
00067 dvector operator*(const dmatrix& m, const dvector& x)
00068  {
00069 #ifdef DIAG
00070      if( heapcheck() == _HEAPCORRUPT )
00071      {
00072         if (ad_printf)
00073           (*ad_printf)( "Entering dvector * dmat dvec is corrupted.\n" );
00074      }
00075      else
00076      {
00077         if (ad_printf)
00078           (*ad_printf)( "Entering dvector * dmat dvec   Heap is OK.\n" );
00079      }
00080 #endif
00081    if (x.indexmin() != m.colmin() || x.indexmax() != m.colmax())
00082    {
00083      cerr << " Incompatible array bounds in "
00084      "dvector  operator * (const dvector& x, const dmatrix& m)\n";
00085      ad_exit(21);
00086    }
00087 
00088    dvector tmp(m.rowmin(),m.rowmax());
00089 
00090    for (int i=m.rowmin(); i<=m.rowmax(); i++)
00091    {
00092      tmp[i]=0;
00093      for (int j=x.indexmin(); j<=x.indexmax(); j++)
00094      {
00095        tmp[i]+=m[i][j]*x[j];
00096      }
00097    }
00098 #ifdef DIAG
00099      if( heapcheck() == _HEAPCORRUPT )
00100      {
00101         if (ad_printf)
00102           (*ad_printf)( "Leaving dvector * dmat dvec is corrupted.\n" );
00103      }
00104      else
00105      {
00106         if (ad_printf)
00107           (*ad_printf)( "Leaving dvector * dmat dvec   Heap is OK.\n" );
00108      }
00109 #endif
00110    return(tmp);
00111  }
00112 
00117  dmatrix  operator * (const dmatrix& m1,const dmatrix& m2 )
00118  {
00119    if (m1.colmin() != m2.rowmin() || m1.colmax() != m2.rowmax())
00120    {
00121      cerr << " Incompatible array bounds in "
00122      "dmatrix  operator * (const dmatrix& x, const dmatrix& m)\n";
00123      ad_exit(21);
00124    }
00125    dmatrix tmp(m1.rowmin(),m1.rowmax(), m2.colmin(), m2.colmax());
00126    double sum;
00127    for (int j=m2.colmin(); j<=m2.colmax(); j++)
00128    {
00129      dvector m2col=column(m2,j);
00130      for (int i=m1.rowmin(); i<=m1.rowmax(); i++)
00131      {
00132        //const dvector& temp_row = m1.elem(i);
00133        sum=m1.elem(i) * m2col;
00134        tmp.elem(i,j)=sum;
00135      }
00136    }
00137    return(tmp);
00138  }
00139 /*
00140 dmatrix operator*(const dmatrix& m1, const dmatrix& m2 )
00141  {
00142    if (m1.colmin() != m2.rowmin() || m1.colmax() != m2.rowmax())
00143    {
00144      cerr << " Incompatible array bounds in "
00145      "dmatrix  operator * (const dmatrix& x, const dmatrix& m)\n";
00146      ad_exit(21);
00147    }
00148 
00149    dmatrix tmp(m1.rowmin(),m1.rowmax(), m2.colmin(), m2.colmax());
00150 
00151    double sum;
00152    const double ** temp_col=
00153          (const double **) malloc(m2.rowsize() * sizeof(double **));
00154    temp_col-=m2.rowmin();
00155 
00156 
00157    for (int j=m2.colmin(); j<=m2.colmax(); j++)
00158    {
00159      for (int k=m2.rowmin(); k<=m2.rowmax(); k++)
00160      {
00161        temp_col[k]=&(m2.elem(k,j));
00162      }
00163 
00164      for (int i=m1.rowmin(); i<=m1.rowmax(); i++)
00165      {
00166        sum=0;
00167 
00168        const dvector& temp_row = m1.elem(i);
00169 
00170        for (k=m1.colmin(); k<=m1.colmax(); k++)
00171        {
00172          sum+=temp_row.elem(k) * *(temp_col[k]);
00173        }
00174        tmp.elem(i,j)=sum;
00175      }
00176    }
00177    temp_col+=m2.rowmin();
00178    free ((char*)temp_col);
00179    return(tmp);
00180  }
00181 */