ADMB Documentation  11.1.2403
 All Classes Files Functions Variables Typedefs Friends Defines
orthply2.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id: orthply2.cpp 1713 2014-03-01 00:26:20Z johnoel $
00003  *
00004  * Author: David Fournier
00005  * Copyright (c) 2008-2012 Regents of the University of California
00006  */
00011 #include <fvar.hpp>
00012 
00017 d4_array orthpoly2(int d1,int d2, int n,int m)
00018 {
00019   d4_array A(0,d1,0,d2,1,n,1,m);
00020   d4_array B(0,d1,0,d2,1,n,1,m);
00021   int alpha,beta,i,j,ii,jj;
00022   int N=(1+d1)*(1+d2);
00023   ivector ialpha(1,N);
00024   ivector ibeta(1,N);
00025 
00026 
00027   for (alpha=0;alpha<=d1;alpha++)
00028   {
00029     for (beta=0;beta<=d2;beta++)
00030     {
00031       for (i=1;i<=n;i++)
00032       {
00033         for (j=1;j<=m;j++)
00034         {
00035 #if defined(USE_DDOUBLE)
00036 #undef double
00037           A(alpha,beta,i,j)=pow(double(i-1)/double(n-1),alpha)*
00038             pow(double(j-1)/double(m-1),beta);
00039 #define double dd_real
00040 #else
00041           A(alpha,beta,i,j)=pow(double(i-1)/double(n-1),alpha)*
00042             pow(double(j-1)/double(m-1),beta);
00043 #endif
00044         }
00045       }
00046     }
00047   }
00048   ii=1;
00049   for (alpha=0;alpha<=d1;alpha++)
00050   {
00051     for (beta=0;beta<=d2;beta++)
00052     {
00053       ialpha(ii)=alpha;
00054       ibeta(ii)=beta;
00055       ii++;
00056     }
00057   }
00058   for (ii=1;ii<=N;ii++)
00059   {
00060     //cout     << "X" << endl;
00061     if (ii>1)
00062     {
00063 //cout << dot(B(ialpha(ii-1),ibeta(ii-1)),A(ialpha(ii),ibeta(ii))) << endl;
00064       //cout << dot(B(ialpha(ii-1),ibeta(ii-1)),
00065        //  A(ialpha(ii),ibeta(ii))/norm(A(ialpha(ii),ibeta(ii)))) << endl;
00066     }
00067     B(ialpha(ii),ibeta(ii))=A(ialpha(ii),ibeta(ii))/
00068       norm(A(ialpha(ii),ibeta(ii)));
00069     //if (ii>1)
00070 // cout << dot(B(ialpha(ii-1),ibeta(ii-1)),B(ialpha(ii),ibeta(ii))) << endl;
00071     //cout     << "Y" << endl;
00072     for (jj=ii+1;jj<=N;jj++)
00073     {
00074       A(ialpha(jj),ibeta(jj))-=
00075         dot(B(ialpha(ii),ibeta(ii)),A(ialpha(jj),ibeta(jj)))*
00076         B(ialpha(ii),ibeta(ii));
00077       //cout << dot(B(ialpha(ii),ibeta(ii)),A(ialpha(jj),ibeta(jj))) << endl;
00078     }
00079     //cout     << "Z" << endl;
00080   }
00081   return B;
00082 }
00083 
00088 double dot(const dmatrix& M,const dmatrix& N)
00089 {
00090   int mmin=M.indexmin();
00091   int mmax=M.indexmax();
00092   if (mmin!=N.indexmin() ||
00093       mmax!=N.indexmax() )
00094   {
00095     cerr << "matrix shapes unequal in"
00096       " double dot(const dmatrix& M,const dmatrix& N)"
00097       << endl;
00098     ad_exit(1);
00099   }
00100   if (M(mmin).indexmin()!=N(mmin).indexmin() ||
00101       M(mmin).indexmax()!=N(mmin).indexmax() )
00102   {
00103     cerr << "matrix shapes unequal in"
00104       " double dot(const dmatrix& M,const dmatrix& N)"
00105       << endl;
00106     ad_exit(1);
00107   }
00108   double ssum=0;
00109   for (int i=mmin;i<=mmax;i++)
00110   {
00111     int jmin=M(i).indexmin();
00112     int jmax=M(i).indexmax();
00113     for (int j=jmin;j<=jmax;j++)
00114     {
00115       ssum+=M(i,j)*N(i,j);
00116     }
00117   }
00118   return ssum;
00119 }