RVL: adjtest.hh Source File

00001 /*************************************************************************
00002 
00003 Copyright Rice University, 2004, 2005, 2006
00004 All rights reserved.
00005 
00006 Permission is hereby granted, free of charge, to any person obtaining a
00007 copy of this software and associated documentation files (the "Software"),
00008 to deal in the Software without restriction, including without limitation
00009 the rights to use, copy, modify, merge, publish, distribute, and/or sell
00010 copies of the Software, and to permit persons to whom the Software is
00011 furnished to do so, provided that the above copyright notice(s) and this
00012 permission notice appear in all copies of the Software and that both the
00013 above copyright notice(s) and this permission notice appear in supporting
00014 documentation.
00015 
00016 THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
00017 IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
00018 FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT OF THIRD PARTY
00019 RIGHTS. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR HOLDERS INCLUDED IN THIS
00020 NOTICE BE LIABLE FOR ANY CLAIM, OR ANY SPECIAL INDIRECT OR CONSEQUENTIAL
00021 DAMAGES, OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR
00022 PROFITS, WHETHER IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS
00023 ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF
00024 THIS SOFTWARE.
00025 
00026 Except as contained in this notice, the name of a copyright holder shall
00027 not be used in advertising or otherwise to promote the sale, use or other
00028 dealings in this Software without prior written authorization of the
00029 copyright holder.
00030 
00031 **************************************************************************/
00032 
00033 #ifndef __RVL_ADJTEST
00034 #define __RVL_ADJTEST
00035 
00036 #include "op.hh"
00037 
00038 namespace RVL {
00039 
00045   template<typename Scalar>
00046   bool AdjointTest(LinearOp<Scalar> const & op,
00047            FunctionObject & randomize,
00048            ostream & str,
00049            int tol=100) {
00050 
00051     typedef typename ScalarFieldTraits<Scalar>::AbsType atype;
00052 
00053     try {
00054 
00055       Vector<Scalar> xin(op.getDomain());
00056       Vector<Scalar> yin(op.getRange());
00057       Vector<Scalar> xout(op.getDomain());
00058       Vector<Scalar> yout(op.getRange());
00059 
00060       xin.eval(randomize);
00061       yin.eval(randomize);
00062   
00063       xout.zero();
00064       yout.zero();
00065   
00066       op.applyOp(xin,yout);
00067       op.applyAdjOp(yin,xout);
00068       Scalar ipdom=xin.inner(xout);
00069       Scalar iprng=yout.inner(yin);
00070       atype denom=yout.norm()*yin.norm();
00071       atype eps = tol*numeric_limits<atype>::epsilon();
00072       // here |Ax| is used as an estimate for |A||x|
00073     
00074       int nd;
00075       if (numeric_precision<atype>()==1) nd = 8;
00076       if (numeric_precision<atype>()==2) nd = 16;
00077       str<<setprecision(nd)<<setiosflags(ios_base::scientific);
00078 
00079       atype aip=abs(ipdom-iprng);
00080       if (aip<eps*denom) {
00081     str<<"AdjointTest"<<endl<<endl;
00082     str<<"adjoint relation holds:";
00083     atype rat;
00084     if (ProtectedDivision<atype>(aip,denom,rat)) {
00085       str<<setprecision(nd)<<setiosflags(ios_base::scientific)<<"<Ax,y>  = "<<iprng<<endl;
00086       str<<setprecision(nd)<<setiosflags(ios_base::scientific)<<"<x,A'y> = "<<ipdom<<endl;
00087       str<<setprecision(nd)<<setiosflags(ios_base::scientific)<<"|Ax||y|  = "<<denom<<endl;
00088       str<<setprecision(nd)<<setiosflags(ios_base::scientific)<<"|<Ax,y>-<x,A'y>| = "<<aip<<endl;
00089       str<<setprecision(nd)<<setiosflags(ios_base::scientific)<<"< 100*macheps*|Ax||y| = "<<eps*denom<<endl;
00090     }
00091     else {
00092       str<<setprecision(nd)<<setiosflags(ios_base::scientific)<<"|<Ax,y>-<x,A'y>|/|Ax||y| = "<<rat<<endl;
00093       str<<setprecision(nd)<<setiosflags(ios::scientific)<<"< 100*macheps = "<<eps<<endl;
00094       str<<"<Ax,y>  = "<<setprecision(nd)<<setiosflags(ios::scientific)<<iprng<<endl;
00095       str<<setprecision(nd)<<setiosflags(ios::scientific)<<"<x,A'y> = "<<ipdom<<endl;
00096       str<<setprecision(nd)<<setiosflags(ios::scientific)<<"|Ax||y|  = "<<denom<<endl;
00097     }
00098     return true;
00099       }
00100       else {
00101     str<<"AdjointTest"<<endl;
00102     str<<"adjoint relation fails: ";
00103     
00104     atype rat;
00105     if (ProtectedDivision<atype>(aip,denom,rat)) {
00106       str<<setprecision(nd)<<setiosflags(ios_base::scientific)<<"<Ax,y>  = "<<iprng<<endl;
00107       str<<setprecision(nd)<<setiosflags(ios_base::scientific)<<"<x,A'y> = "<<ipdom<<endl;
00108       str<<setprecision(nd)<<setiosflags(ios_base::scientific)<<"|Ax||y|  = "<<denom<<endl;
00109       str<<setprecision(nd)<<setiosflags(ios_base::scientific)<<"|<Ax,y>-<x,A'y>| = "<<aip<<endl;
00110       str<<setprecision(nd)<<setiosflags(ios_base::scientific)<<">= 100*macheps*|Ax||y| = "<<eps*denom<<endl;     
00111     }
00112     else {
00113       str<<setprecision(nd)<<setiosflags(ios_base::scientific)<<"|<Ax,y>-<x,A'y>|/|Ax||y| = "<<rat<<endl;
00114       str<<setprecision(nd)<<setiosflags(ios_base::scientific)<<">= 100*macheps = "<<eps<<endl;     
00115       str<<setprecision(nd)<<setiosflags(ios_base::scientific)<<"<Ax,y>  = "<<iprng<<endl;
00116       str<<setprecision(nd)<<setiosflags(ios_base::scientific)<<"<x,A'y> = "<<ipdom<<endl;
00117       str<<setprecision(nd)<<setiosflags(ios_base::scientific)<<"|Ax||y|  = "<<denom<<endl;
00118     }
00119       }
00120       return false;
00121     }
00122     catch (RVLException & e) {
00123       e<<"\ncalled from AdjointTest\n";
00124       throw e;
00125     }
00126   }
00127 }
00128 
00129 #endif