ADMB Documentation  11.1.2192
 All Classes Files Functions Variables Typedefs Friends Defines
dmat1.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id: dmat1.cpp 1711 2014-02-28 22:46:44Z 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        sum=0;
00133        //const dvector& temp_row = m1.elem(i);
00134        sum=m1.elem(i) * m2col;
00135        tmp.elem(i,j)=sum;
00136      }
00137    }
00138    return(tmp);
00139  }
00140 /*
00141 dmatrix operator*(const dmatrix& m1, const dmatrix& m2 )
00142  {
00143    if (m1.colmin() != m2.rowmin() || m1.colmax() != m2.rowmax())
00144    {
00145      cerr << " Incompatible array bounds in "
00146      "dmatrix  operator * (const dmatrix& x, const dmatrix& m)\n";
00147      ad_exit(21);
00148    }
00149 
00150    dmatrix tmp(m1.rowmin(),m1.rowmax(), m2.colmin(), m2.colmax());
00151 
00152    double sum;
00153    const double ** temp_col=
00154          (const double **) malloc(m2.rowsize() * sizeof(double **));
00155    temp_col-=m2.rowmin();
00156 
00157 
00158    for (int j=m2.colmin(); j<=m2.colmax(); j++)
00159    {
00160      for (int k=m2.rowmin(); k<=m2.rowmax(); k++)
00161      {
00162        temp_col[k]=&(m2.elem(k,j));
00163      }
00164 
00165      for (int i=m1.rowmin(); i<=m1.rowmax(); i++)
00166      {
00167        sum=0;
00168 
00169        const dvector& temp_row = m1.elem(i);
00170 
00171        for (k=m1.colmin(); k<=m1.colmax(); k++)
00172        {
00173          sum+=temp_row.elem(k) * *(temp_col[k]);
00174        }
00175        tmp.elem(i,j)=sum;
00176      }
00177    }
00178    temp_col+=m2.rowmin();
00179    free ((char*)temp_col);
00180    return(tmp);
00181  }
00182 */