ADMB Documentation  11.1.1903
 All Classes Files Functions Variables Typedefs Friends Defines
ad_atlas.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id: ad_atlas.cpp 1136 2013-08-05 19:40:39Z johnoel $
00003  *
00004  * Author: David Fournier
00005  * Copyright (c) 2008-2012 Regents of the University of California
00006  */
00011 #include <admodel.h>
00012 
00013 #if defined USE_ATLAS
00014 
00015 #include <clapack.h>
00016 
00021 dvector atlas_solve_spd(const dmatrix & M, const dvector & x)
00022 {
00023   int mmin=M.indexmin();
00024   int mmax=M.indexmax();
00025   if (mmin != M(mmin).indexmin() ||
00026      mmax != M(mmin).indexmax())
00027   {
00028     cerr << "Matrix not square in dvector atlas_solve_spd"
00029          << endl;
00030     ad_exit(1);
00031   }
00032   if (mmin != x.indexmin() || mmax != x.indexmax())
00033   {
00034     cerr << "Incompatible matrix and vector sizes in dvector atlas_solve_spd"
00035          << endl;
00036     ad_exit(1);
00037   }
00038   dvector v(mmin,mmax);
00039   int sz=mmax-mmin+1;
00040   dvector M1(1,sz*sz);
00041   int ii=1;
00042   int i,j;
00043   for (i=mmin;i<=mmax;i++)
00044   {
00045     for (j=mmin;j<=mmax;j++)
00046     {
00047       M1(ii++)=M(i,j);
00048     }
00049   }
00050   v=x;
00051   double *Ap= &(M1(1));
00052   double *X = &(v(mmin)); const int incX=1;
00053 
00054   const enum CBLAS_ORDER Order=CblasRowMajor;
00055   const enum CBLAS_UPLO Uplo=CblasLower;
00056   const enum CBLAS_TRANSPOSE TransA=CblasNoTrans;
00057   const enum CBLAS_DIAG Diag=CblasNonUnit;
00058 
00059   int retr=clapack_dposv(Order, Uplo, sz,1, Ap, sz, X, sz);
00060   return v;
00061 }
00062 
00067 dvector atlas_solve_spd(const dmatrix & M, const dvector & x, int& ierr)
00068 {
00069   int mmin=M.indexmin();
00070   int mmax=M.indexmax();
00071   if (mmin != M(mmin).indexmin() ||
00072      mmax != M(mmin).indexmax())
00073   {
00074     cerr << "Matrix not square in dvector atlas_solve_spd"
00075          << endl;
00076     ad_exit(1);
00077   }
00078   if (mmin != x.indexmin() || mmax != x.indexmax())
00079   {
00080     cerr << "Incompatible matrix and vector sizes in dvector atlas_solve_spd"
00081          << endl;
00082     ad_exit(1);
00083   }
00084   dvector v(mmin,mmax);
00085   int sz=mmax-mmin+1;
00086   dvector M1(1,sz*sz);
00087   int ii=1;
00088   int i,j;
00089   for (i=mmin;i<=mmax;i++)
00090   {
00091     for (j=mmin;j<=mmax;j++)
00092     {
00093       M1(ii++)=M(i,j);
00094     }
00095   }
00096   v=x;
00097   double *Ap= &(M1(1));
00098   double *X = &(v(mmin)); const int incX=1;
00099 
00100   const enum CBLAS_ORDER Order=CblasRowMajor;
00101   const enum CBLAS_UPLO Uplo=CblasLower;
00102   const enum CBLAS_TRANSPOSE TransA=CblasNoTrans;
00103   const enum CBLAS_DIAG Diag=CblasNonUnit;
00104 
00105   int retr=clapack_dposv(Order, Uplo, sz,1, Ap, sz, X, sz);
00106   ierr=retr;
00107   return v;
00108 }
00109 
00114 dmatrix atlas_solve_spd(const dmatrix & M, const dmatrix & N)
00115 {
00116   int mmin=M.indexmin();
00117   int mmax=M.indexmax();
00118   int nmin=N.indexmin();
00119   int nmax=N.indexmax();
00120   if (mmin != M(mmin).indexmin() ||
00121      mmax != M(mmin).indexmax())
00122   {
00123     cerr << "Matrix not square in dvector atlas_solve_spd"
00124          << endl;
00125     ad_exit(1);
00126   }
00127   if (mmin != N(nmin).indexmin() || mmax != N(nmin).indexmax() )
00128   {
00129     cerr << "Incompatible matrix and vector sizes in dmatrix atlas_solve_spd"
00130          << endl;
00131     ad_exit(1);
00132   }
00133   int sz=mmax-mmin+1;
00134   int szn=nmax-nmin+1;
00135   dvector M1(1,sz*sz);
00136   dvector N1(1,sz*szn);
00137   int ii=1;
00138   int i,j;
00139   for (i=mmin;i<=mmax;i++)
00140   {
00141     for (j=mmin;j<=mmax;j++)
00142     {
00143       M1(ii++)=M(i,j);
00144     }
00145   }
00146   ii=1;
00147   for (i=nmin;i<=nmax;i++)
00148   {
00149     for (j=mmin;j<=mmax;j++)
00150     {
00151       N1(ii++)=N(i,j);
00152     }
00153   }
00154   double *Ap= &(M1(1));
00155   double *X = &(N1(1)); const int incX=1;
00156 
00157   const enum CBLAS_ORDER Order=CblasRowMajor;
00158   const enum CBLAS_UPLO Uplo=CblasLower;
00159   const enum CBLAS_TRANSPOSE TransA=CblasNoTrans;
00160   const enum CBLAS_DIAG Diag=CblasNonUnit;
00161 
00162   //int retr=clapack_dposv(Order, Uplo, sz,1  , Ap, sz, X, sz);
00163   int retr=clapack_dposv(Order, Uplo, sz,szn, Ap, sz, X, sz);
00164   dmatrix tmp(nmin,nmax,mmin,mmax);
00165 
00166   ii=0;
00167   for (i=nmin;i<=nmax;i++)
00168   {
00169     for (j=mmin;j<=mmax;j++)
00170     {
00171       tmp(i,j)=X[ii++];
00172     }
00173   }
00174   return tmp;
00175 }
00176 
00181 dmatrix atlas_solve_spd_trans(const dmatrix & M, const dmatrix & N)
00182 {
00183   int mmin=M.indexmin();
00184   int mmax=M.indexmax();
00185   if (mmin != M(mmin).indexmin() ||
00186      mmax != M(mmin).indexmax())
00187   {
00188     cerr << "Matrix not square in dvector atlas_solve_spd"
00189          << endl;
00190     ad_exit(1);
00191   }
00192   if (mmin != N.indexmin() || mmax != N.indexmax() )
00193   {
00194     cerr << "Incompatible matrix and vector sizes in dmatrix atlas_solve_spd"
00195          << endl;
00196     ad_exit(1);
00197   }
00198   int nmin=N(mmin).indexmin();
00199   int nmax=N(mmin).indexmax();
00200   int sz=mmax-mmin+1;
00201   int szn=nmax-nmin+1;
00202   dvector M1(1,sz*sz);
00203   dvector N1(1,sz*szn);
00204   int ii=1;
00205   int i,j;
00206   for (i=mmin;i<=mmax;i++)
00207   {
00208     for (j=mmin;j<=mmax;j++)
00209     {
00210       M1(ii++)=M(i,j);
00211     }
00212   }
00213   ii=1;
00214   for (j=nmin;j<=nmax;j++)
00215   {
00216     for (i=mmin;i<=mmax;i++)
00217     {
00218       N1(ii++)=N(i,j);
00219     }
00220   }
00221   double *Ap= &(M1(1));
00222   double *X = &(N1(1)); const int incX=1;
00223 
00224   const enum CBLAS_ORDER Order=CblasRowMajor;
00225   const enum CBLAS_UPLO Uplo=CblasLower;
00226   const enum CBLAS_TRANSPOSE TransA=CblasNoTrans;
00227   const enum CBLAS_DIAG Diag=CblasNonUnit;
00228 
00229   //int retr=clapack_dposv(Order, Uplo, sz,1  , Ap, sz, X, sz);
00230   int retr=clapack_dposv(Order, Uplo, sz,szn, Ap, sz, X, sz);
00231   dmatrix tmp(mmin,mmax,nmin,nmax);
00232 
00233   ii=0;
00234   for (j=nmin;j<=nmax;j++)
00235   {
00236     for (i=mmin;i<=mmax;i++)
00237     {
00238       tmp(i,j)=X[ii++];
00239     }
00240   }
00241   return tmp;
00242 }
00243 #endif //#if defined USE_ATLAS