<html><head><meta http-equiv="Content-Type" content="text/html;charset=iso-8859-1">
<title>HippoDraw Class Library</title>
<link href="doxygen.css" rel="stylesheet" type="text/css">
<link href="tabs.css" rel="stylesheet" type="text/css">
<link rel="shortcut icon" href="hippoApp.ico">
</head><body>
<!-- Generated by Doxygen 1.6.2 -->
<div class="navigation" id="top">
<div class="tabs">
<ul>
<li><a href="index.html"><span>Main Page</span></a></li>
<li><a href="pages.html"><span>Related Pages</span></a></li>
<li><a href="namespaces.html"><span>Namespaces</span></a></li>
<li><a href="annotated.html"><span>Classes</span></a></li>
<li class="current"><a href="files.html"><span>Files</span></a></li>
<li><a href="dirs.html"><span>Directories</span></a></li>
</ul>
</div>
<div class="tabs">
<ul>
<li><a href="files.html"><span>File List</span></a></li>
<li><a href="globals.html"><span>File Members</span></a></li>
</ul>
</div>
<div class="navpath"><a class="el" href="dir_6633238f5107fd0d748d9093e08c459d.html">minimizers</a>
</div>
</div>
<div class="contents">
<h1>NumLinAlg.cxx</h1><a href="NumLinAlg_8cxx.html">Go to the documentation of this file.</a><div class="fragment"><pre class="fragment"><a name="l00001"></a>00001
<a name="l00014"></a>00014 <span class="preprocessor">#include "<a class="code" href="NumLinAlg_8h.html" title="Collection of linear algebra functions.">NumLinAlg.h</a>"</span>
<a name="l00015"></a>00015
<a name="l00016"></a>00016 <span class="preprocessor">#include <fstream></span>
<a name="l00017"></a>00017 <span class="preprocessor">#include <iostream></span>
<a name="l00018"></a>00018 <span class="preprocessor">#include <iomanip></span>
<a name="l00019"></a>00019 <span class="preprocessor">#include <cmath></span>
<a name="l00020"></a>00020 <span class="preprocessor">#include <cfloat></span>
<a name="l00021"></a>00021 <span class="preprocessor">#include <cassert></span>
<a name="l00022"></a>00022 <span class="preprocessor">#include <cstdlib></span>
<a name="l00023"></a>00023
<a name="l00024"></a>00024 <span class="keyword">using</span> std::ofstream;
<a name="l00025"></a>00025 <span class="keyword">using</span> std::ifstream;
<a name="l00026"></a>00026 <span class="keyword">using</span> std::setw;
<a name="l00027"></a>00027 <span class="keyword">using</span> std::setprecision;
<a name="l00028"></a>00028 <span class="keyword">using</span> std::cout;
<a name="l00029"></a>00029 <span class="keyword">using</span> std::endl;
<a name="l00030"></a>00030 <span class="keyword">using</span> std::vector;
<a name="l00031"></a>00031 <span class="keyword">using</span> std::abs;
<a name="l00032"></a>00032
<a name="l00033"></a>00033 <span class="keyword">namespace </span>hippodraw {
<a name="l00034"></a>00034 <span class="keyword">namespace </span>Numeric {
<a name="l00035"></a>00035
<a name="l00036"></a>00036 std::vector< double >
<a name="l00037"></a>00037 <a class="code" href="namespacehippodraw_1_1Numeric.html#aab0ab86ece0308067f6955d5fbd36349" title="Given two vectors x and y this function performs operation z = x + y.">operator + </a>( <span class="keyword">const</span> std::vector< double >& x,
<a name="l00038"></a>00038 <span class="keyword">const</span> std::vector< double >& y )
<a name="l00039"></a>00039 {
<a name="l00040"></a>00040 <span class="keywordtype">int</span> nrows = x.size();
<a name="l00041"></a>00041 vector< double > z;
<a name="l00042"></a>00042
<a name="l00043"></a>00043 z.resize( nrows, 0.0 );
<a name="l00044"></a>00044
<a name="l00045"></a>00045 <span class="keywordflow">for</span> ( <span class="keywordtype">int</span> i = 0; i < nrows; i++)
<a name="l00046"></a>00046 z[i] = x[i] + y[i];
<a name="l00047"></a>00047
<a name="l00048"></a>00048 <span class="keywordflow">return</span> z;
<a name="l00049"></a>00049 }
<a name="l00050"></a>00050
<a name="l00051"></a>00051 std::vector< double >
<a name="l00052"></a>00052 <a class="code" href="namespacehippodraw_1_1Numeric.html#a83635d69c8d1388f6f5857be75105d43" title="Given two vectors x and y this function performs operation z = x - y.">operator - </a>( <span class="keyword">const</span> std::vector< double >& x,
<a name="l00053"></a>00053 <span class="keyword">const</span> std::vector< double >& y )
<a name="l00054"></a>00054 {
<a name="l00055"></a>00055 <span class="keywordtype">int</span> nrows = x.size();
<a name="l00056"></a>00056 vector< double > z;
<a name="l00057"></a>00057
<a name="l00058"></a>00058 z.resize( nrows, 0.0 );
<a name="l00059"></a>00059
<a name="l00060"></a>00060 <span class="keywordflow">for</span> ( <span class="keywordtype">int</span> i = 0; i < nrows; i++)
<a name="l00061"></a>00061 z[i] = x[i] - y[i];
<a name="l00062"></a>00062
<a name="l00063"></a>00063 <span class="keywordflow">return</span> z;
<a name="l00064"></a>00064 }
<a name="l00065"></a>00065
<a name="l00066"></a>00066 std::vector< double >
<a name="l00067"></a>00067 <a class="code" href="namespacehippodraw_1_1Numeric.html#a83635d69c8d1388f6f5857be75105d43" title="Given two vectors x and y this function performs operation z = x - y.">operator - </a>( <span class="keyword">const</span> std::vector< double >& y )
<a name="l00068"></a>00068 {
<a name="l00069"></a>00069 <span class="keywordtype">int</span> nrows = y.size();
<a name="l00070"></a>00070 vector< double > z;
<a name="l00071"></a>00071
<a name="l00072"></a>00072 z.resize( nrows, 0.0 );
<a name="l00073"></a>00073
<a name="l00074"></a>00074 <span class="keywordflow">for</span> ( <span class="keywordtype">int</span> i = 0; i < nrows; i++ )
<a name="l00075"></a>00075 z[i] = - y[i];
<a name="l00076"></a>00076
<a name="l00077"></a>00077 <span class="keywordflow">return</span> z;
<a name="l00078"></a>00078 }
<a name="l00079"></a>00079
<a name="l00080"></a>00080 std::vector< vector< double > >
<a name="l00081"></a>00081 <a class="code" href="namespacehippodraw_1_1Numeric.html#aab0ab86ece0308067f6955d5fbd36349" title="Given two vectors x and y this function performs operation z = x + y.">operator + </a>( <span class="keyword">const</span> std::vector< std::vector< double > >&A,
<a name="l00082"></a>00082 <span class="keyword">const</span> std::vector< std::vector< double > >&B )
<a name="l00083"></a>00083 {
<a name="l00084"></a>00084 <span class="keywordtype">int</span> nrows = A.size();
<a name="l00085"></a>00085 <span class="keywordtype">int</span> ncols = A[0].size();
<a name="l00086"></a>00086
<a name="l00087"></a>00087 vector< vector< double > > C( nrows );
<a name="l00088"></a>00088 <span class="keywordflow">for</span>( <span class="keywordtype">int</span> i = 0; i < nrows; i++ )
<a name="l00089"></a>00089 C[i].resize( ncols, 0.0 );
<a name="l00090"></a>00090
<a name="l00091"></a>00091 <span class="keywordflow">for</span> (<span class="keywordtype">int</span> i = 0; i < nrows; i++)
<a name="l00092"></a>00092 <span class="keywordflow">for</span> (<span class="keywordtype">int</span> j = 0; j < ncols; j++)
<a name="l00093"></a>00093 C[i][j] = A[i][j] + B[i][j];
<a name="l00094"></a>00094
<a name="l00095"></a>00095 <span class="keywordflow">return</span> C;
<a name="l00096"></a>00096 }
<a name="l00097"></a>00097
<a name="l00098"></a>00098 std::vector< vector< double > >
<a name="l00099"></a>00099 <a class="code" href="namespacehippodraw_1_1Numeric.html#a83635d69c8d1388f6f5857be75105d43" title="Given two vectors x and y this function performs operation z = x - y.">operator - </a>( <span class="keyword">const</span> std::vector< std::vector< double > >&A,
<a name="l00100"></a>00100 <span class="keyword">const</span> std::vector< std::vector< double > >&B )
<a name="l00101"></a>00101 {
<a name="l00102"></a>00102 <span class="keywordtype">int</span> nrows = A.size();
<a name="l00103"></a>00103 <span class="keywordtype">int</span> ncols = A[0].size();
<a name="l00104"></a>00104
<a name="l00105"></a>00105 vector< vector< double > > C( nrows );
<a name="l00106"></a>00106 <span class="keywordflow">for</span>( <span class="keywordtype">int</span> i = 0; i < nrows; i++ )
<a name="l00107"></a>00107 C[i].resize( ncols, 0.0 );
<a name="l00108"></a>00108
<a name="l00109"></a>00109 <span class="keywordflow">for</span> (<span class="keywordtype">int</span> i = 0; i < nrows; i++)
<a name="l00110"></a>00110 <span class="keywordflow">for</span> (<span class="keywordtype">int</span> j = 0; j < ncols; j++)
<a name="l00111"></a>00111 C[i][j] = A[i][j] - B[i][j];
<a name="l00112"></a>00112
<a name="l00113"></a>00113 <span class="keywordflow">return</span> C;
<a name="l00114"></a>00114 }
<a name="l00115"></a>00115
<a name="l00116"></a>00116
<a name="l00117"></a>00117 std::vector< double >
<a name="l00118"></a>00118 <a class="code" href="namespacehippodraw_1_1Numeric.html#a25db1c1e5a7b83f6b4b0ea16793808fd" title="Given a scalar a and a vector x this function performs operation y = ax.">operator * </a>( <span class="keywordtype">double</span> a, <span class="keyword">const</span> std::vector< double >& x )
<a name="l00119"></a>00119 {
<a name="l00120"></a>00120 <span class="keywordtype">int</span> nrows = x.size();
<a name="l00121"></a>00121 vector< double > y;
<a name="l00122"></a>00122
<a name="l00123"></a>00123 y.resize( nrows, 0.0 );
<a name="l00124"></a>00124
<a name="l00125"></a>00125 <span class="keywordflow">for</span> ( <span class="keywordtype">int</span> i = 0; i < nrows; i++)
<a name="l00126"></a>00126 y[i] = a * x[i];
<a name="l00127"></a>00127
<a name="l00128"></a>00128 <span class="keywordflow">return</span> y;
<a name="l00129"></a>00129 }
<a name="l00130"></a>00130
<a name="l00131"></a>00131 std::vector< double >
<a name="l00132"></a>00132 <a class="code" href="namespacehippodraw_1_1Numeric.html#a1f5254a9fe7eb0973dd1eff892ac39a0" title="Given a scalar and a vector x this function performs operation y = x/a.">operator / </a>( <span class="keyword">const</span> std::vector< double >& x, <span class="keywordtype">double</span> a )
<a name="l00133"></a>00133 {
<a name="l00134"></a>00134 <span class="keywordtype">int</span> nrows = x.size();
<a name="l00135"></a>00135 vector< double > y;
<a name="l00136"></a>00136
<a name="l00137"></a>00137 assert( abs( a ) > DBL_EPSILON );
<a name="l00138"></a>00138
<a name="l00139"></a>00139 y.resize( nrows, 0.0 );
<a name="l00140"></a>00140
<a name="l00141"></a>00141 <span class="keywordflow">for</span> ( <span class="keywordtype">int</span> i = 0; i < nrows; i++)
<a name="l00142"></a>00142 y[i] = x[i] / a;
<a name="l00143"></a>00143
<a name="l00144"></a>00144 <span class="keywordflow">return</span> y;
<a name="l00145"></a>00145 }
<a name="l00146"></a>00146
<a name="l00147"></a>00147 std::vector< std::vector< double > >
<a name="l00148"></a>00148 <a class="code" href="namespacehippodraw_1_1Numeric.html#a25db1c1e5a7b83f6b4b0ea16793808fd" title="Given a scalar a and a vector x this function performs operation y = ax.">operator * </a>( <span class="keywordtype">double</span> a, <span class="keyword">const</span> std::vector< std::vector< double > >&A )
<a name="l00149"></a>00149 {
<a name="l00150"></a>00150 <span class="keywordtype">int</span> nrows = A.size();
<a name="l00151"></a>00151 <span class="keywordtype">int</span> ncols = A[0].size();
<a name="l00152"></a>00152
<a name="l00153"></a>00153 vector< vector< double > > B( nrows );
<a name="l00154"></a>00154 <span class="keywordflow">for</span>( <span class="keywordtype">int</span> i = 0; i < nrows; i++ )
<a name="l00155"></a>00155 B[i].resize( ncols, 0.0 );
<a name="l00156"></a>00156
<a name="l00157"></a>00157 <span class="keywordflow">for</span> (<span class="keywordtype">int</span> i = 0; i < nrows; i++)
<a name="l00158"></a>00158 <span class="keywordflow">for</span> (<span class="keywordtype">int</span> j = 0; j < ncols; j++)
<a name="l00159"></a>00159 B[i][j] = a * A[i][j];
<a name="l00160"></a>00160
<a name="l00161"></a>00161 <span class="keywordflow">return</span> B;
<a name="l00162"></a>00162 }
<a name="l00163"></a>00163
<a name="l00164"></a>00164 std::vector< std::vector< double > >
<a name="l00165"></a>00165 <a class="code" href="namespacehippodraw_1_1Numeric.html#a1f5254a9fe7eb0973dd1eff892ac39a0" title="Given a scalar and a vector x this function performs operation y = x/a.">operator / </a>( <span class="keyword">const</span> std::vector< std::vector< double > >& A, <span class="keywordtype">double</span> a )
<a name="l00166"></a>00166 {
<a name="l00167"></a>00167 <span class="keywordtype">int</span> nrows = A.size();
<a name="l00168"></a>00168 <span class="keywordtype">int</span> ncols = A[0].size();
<a name="l00169"></a>00169
<a name="l00170"></a>00170 assert( abs( a ) > DBL_EPSILON );
<a name="l00171"></a>00171
<a name="l00172"></a>00172 vector< vector< double > > B( nrows );
<a name="l00173"></a>00173 <span class="keywordflow">for</span>( <span class="keywordtype">int</span> i = 0; i < nrows; i++ )
<a name="l00174"></a>00174 B[i].resize( ncols, 0.0 );
<a name="l00175"></a>00175
<a name="l00176"></a>00176 <span class="keywordflow">for</span> (<span class="keywordtype">int</span> i = 0; i < nrows; i++)
<a name="l00177"></a>00177 <span class="keywordflow">for</span> (<span class="keywordtype">int</span> j = 0; j < ncols; j++)
<a name="l00178"></a>00178 B[i][j] = A[i][j]/a;
<a name="l00179"></a>00179
<a name="l00180"></a>00180 <span class="keywordflow">return</span> B;
<a name="l00181"></a>00181 }
<a name="l00182"></a>00182
<a name="l00183"></a>00183 std::vector< double >
<a name="l00184"></a>00184 <a class="code" href="namespacehippodraw_1_1Numeric.html#a25db1c1e5a7b83f6b4b0ea16793808fd" title="Given a scalar a and a vector x this function performs operation y = ax.">operator * </a>( <span class="keyword">const</span> std::vector< std::vector< double > >& A,
<a name="l00185"></a>00185 <span class="keyword">const</span> std::vector< double >& x )
<a name="l00186"></a>00186 {
<a name="l00187"></a>00187 <span class="keywordtype">int</span> nrows = A.size();
<a name="l00188"></a>00188 <span class="keywordtype">int</span> ncols = A[0].size();
<a name="l00189"></a>00189 vector< double > y;
<a name="l00190"></a>00190
<a name="l00191"></a>00191 y.resize( nrows, 0.0 );
<a name="l00192"></a>00192
<a name="l00193"></a>00193 <span class="keywordflow">for</span> (<span class="keywordtype">int</span> i = 0; i < nrows; i++)
<a name="l00194"></a>00194 <span class="keywordflow">for</span> (<span class="keywordtype">int</span> j = 0; j < ncols; j++)
<a name="l00195"></a>00195 y[i] += A[i][j] * x[j];
<a name="l00196"></a>00196
<a name="l00197"></a>00197 <span class="keywordflow">return</span> y;
<a name="l00198"></a>00198 }
<a name="l00199"></a>00199
<a name="l00200"></a>00200 std::vector< double >
<a name="l00201"></a>00201 <a class="code" href="namespacehippodraw_1_1Numeric.html#a25db1c1e5a7b83f6b4b0ea16793808fd" title="Given a scalar a and a vector x this function performs operation y = ax.">operator * </a>( <span class="keyword">const</span> std::vector< double >& x,
<a name="l00202"></a>00202 <span class="keyword">const</span> std::vector< std::vector< double > >& A )
<a name="l00203"></a>00203 {
<a name="l00204"></a>00204 <span class="keywordtype">int</span> nrows = A.size();
<a name="l00205"></a>00205 <span class="keywordtype">int</span> ncols = A[0].size();
<a name="l00206"></a>00206 vector< double > y;
<a name="l00207"></a>00207
<a name="l00208"></a>00208 y.resize( ncols, 0.0 );
<a name="l00209"></a>00209
<a name="l00210"></a>00210 <span class="keywordflow">for</span> (<span class="keywordtype">int</span> j = 0; j < ncols; j++)
<a name="l00211"></a>00211 <span class="keywordflow">for</span> (<span class="keywordtype">int</span> i = 0; i < nrows; i++)
<a name="l00212"></a>00212 y[j] += A[i][j] * x[i];
<a name="l00213"></a>00213
<a name="l00214"></a>00214 <span class="keywordflow">return</span> y;
<a name="l00215"></a>00215 }
<a name="l00216"></a>00216
<a name="l00217"></a>00217 std::vector< vector< double > >
<a name="l00218"></a>00218 <a class="code" href="namespacehippodraw_1_1Numeric.html#a25db1c1e5a7b83f6b4b0ea16793808fd" title="Given a scalar a and a vector x this function performs operation y = ax.">operator * </a>( <span class="keyword">const</span> std::vector< std::vector< double > >&A,
<a name="l00219"></a>00219 <span class="keyword">const</span> std::vector< std::vector< double > >&B )
<a name="l00220"></a>00220 {
<a name="l00221"></a>00221 <span class="keywordtype">int</span> m = A.size();
<a name="l00222"></a>00222 <span class="keywordtype">int</span> r = A[0].size();
<a name="l00223"></a>00223 <span class="keywordtype">int</span> n = B[0].size();
<a name="l00224"></a>00224
<a name="l00225"></a>00225 vector< vector< double > > C( m );
<a name="l00226"></a>00226 <span class="keywordflow">for</span>( <span class="keywordtype">int</span> i = 0; i < m; i++ )
<a name="l00227"></a>00227 C[i].resize( n, 0.0 );
<a name="l00228"></a>00228
<a name="l00229"></a>00229 <span class="keywordflow">for</span> (<span class="keywordtype">int</span> i = 0; i < m; i++)
<a name="l00230"></a>00230 <span class="keywordflow">for</span> (<span class="keywordtype">int</span> j = 0; j < n; j++)
<a name="l00231"></a>00231 <span class="keywordflow">for</span> (<span class="keywordtype">int</span> k = 0; k < r; k++)
<a name="l00232"></a>00232 C[i][j] += A[i][k] * B[k][j];
<a name="l00233"></a>00233
<a name="l00234"></a>00234 <span class="keywordflow">return</span> C;
<a name="l00235"></a>00235 }
<a name="l00236"></a>00236
<a name="l00237"></a>00237 <span class="keywordtype">double</span>
<a name="l00238"></a>00238 <a class="code" href="namespacehippodraw_1_1Numeric.html#a329ce3200d86ee77e2b9721db28388ee" title="Computes the dot or the inner product of two vectors(i.e.">innerProduct</a>( <span class="keyword">const</span> std::vector< double > & a,
<a name="l00239"></a>00239 <span class="keyword">const</span> std::vector< double > & b )
<a name="l00240"></a>00240 {
<a name="l00241"></a>00241 <span class="keywordtype">double</span> prod = 0.0;
<a name="l00242"></a>00242
<a name="l00243"></a>00243 <span class="keywordflow">for</span> ( <span class="keywordtype">unsigned</span> <span class="keywordtype">int</span> i = 0; i < a.size(); i++ )
<a name="l00244"></a>00244 prod += a[i] * b[i];
<a name="l00245"></a>00245
<a name="l00246"></a>00246 <span class="keywordflow">return</span> prod;
<a name="l00247"></a>00247 }
<a name="l00248"></a>00248
<a name="l00249"></a>00249
<a name="l00250"></a>00250 vector< vector< double > >
<a name="l00251"></a>00251 <a class="code" href="namespacehippodraw_1_1Numeric.html#af9cb017e1c9362831eccadf6bceac017" title="Computes the outer product of two vectors (i.e.">outerProduct</a> ( <span class="keyword">const</span> std::vector< double >& a,
<a name="l00252"></a>00252 <span class="keyword">const</span> std::vector< double >& b )
<a name="l00253"></a>00253 {
<a name="l00254"></a>00254 vector< vector< double> > C( a.size() );
<a name="l00255"></a>00255 <span class="keywordflow">for</span>( <span class="keywordtype">unsigned</span> <span class="keywordtype">int</span> i = 0; i < a.size(); i++ )
<a name="l00256"></a>00256 C[i].resize( b.size() );
<a name="l00257"></a>00257
<a name="l00258"></a>00258 <span class="keywordflow">for</span>( <span class="keywordtype">unsigned</span> <span class="keywordtype">int</span> i = 0; i < a.size(); i++ )
<a name="l00259"></a>00259 <span class="keywordflow">for</span>( <span class="keywordtype">unsigned</span> <span class="keywordtype">int</span> j = 0; j < b.size(); j++ )
<a name="l00260"></a>00260 C[i][j] = a[i] * b[j];
<a name="l00261"></a>00261
<a name="l00262"></a>00262 <span class="keywordflow">return</span> C;
<a name="l00263"></a>00263 }
<a name="l00264"></a>00264
<a name="l00265"></a>00265
<a name="l00266"></a>00266 <span class="keywordtype">double</span>
<a name="l00267"></a>00267 <a class="code" href="namespacehippodraw_1_1Numeric.html#ab6ec8f5900d139d769c7e4aaa8e63f8a" title="Calculates the vector-matrix-vector product x&#39;*A*x.">quadraticProduct</a>( <span class="keyword">const</span> std::vector< std::vector< double > > & A,
<a name="l00268"></a>00268 <span class="keyword">const</span> std::vector< double > x )
<a name="l00269"></a>00269 {
<a name="l00270"></a>00270 <span class="keywordtype">double</span> prod = 0.0;
<a name="l00271"></a>00271 <span class="keywordtype">unsigned</span> <span class="keywordtype">int</span> n = A.size();
<a name="l00272"></a>00272
<a name="l00273"></a>00273 assert ( x.size() == n );
<a name="l00274"></a>00274
<a name="l00275"></a>00275 <span class="keywordflow">for</span> ( <span class="keywordtype">unsigned</span> <span class="keywordtype">int</span> i = 0; i < n; i++ )
<a name="l00276"></a>00276 <span class="keywordflow">for</span> ( <span class="keywordtype">unsigned</span> <span class="keywordtype">int</span> j = 0; j < n; j++ )
<a name="l00277"></a>00277 prod += x[i] * A[i][j] * x[j];
<a name="l00278"></a>00278
<a name="l00279"></a>00279 <span class="keywordflow">return</span> prod;
<a name="l00280"></a>00280 }
<a name="l00281"></a>00281
<a name="l00282"></a>00282
<a name="l00283"></a>00283 <span class="keywordtype">double</span>
<a name="l00284"></a>00284 <a class="code" href="namespacehippodraw_1_1Numeric.html#a2fc63b8199da4a6f45b4a9590fbab38f" title="Computes the two norm of the vector.">norm</a> ( <span class="keyword">const</span> std::vector< double > & a )
<a name="l00285"></a>00285 {
<a name="l00286"></a>00286 <span class="keywordflow">return</span> sqrt( <a class="code" href="namespacehippodraw_1_1Numeric.html#a329ce3200d86ee77e2b9721db28388ee" title="Computes the dot or the inner product of two vectors(i.e.">innerProduct</a>( a, a ) );
<a name="l00287"></a>00287 }
<a name="l00288"></a>00288
<a name="l00289"></a>00289
<a name="l00290"></a>00290 <span class="keywordtype">int</span>
<a name="l00291"></a>00291 <a class="code" href="namespacehippodraw_1_1Numeric.html#aecbdf19ba7109e37b7e89e490aaa453a" title="The subroutine which does cholesky factorization of a given Symmetric positive definite...">cholFactor</a> ( std::vector < std::vector< double > > & A )
<a name="l00292"></a>00292 {
<a name="l00293"></a>00293 <span class="keywordtype">unsigned</span> <span class="keywordtype">int</span> n = A.size();
<a name="l00294"></a>00294
<a name="l00295"></a>00295 <span class="keywordflow">for</span> ( <span class="keywordtype">unsigned</span> <span class="keywordtype">int</span> i = 0; i < n ; ++i )
<a name="l00296"></a>00296 <span class="keywordflow">for</span> ( <span class="keywordtype">unsigned</span> <span class="keywordtype">int</span> j = 0; j <= i ; ++j)
<a name="l00297"></a>00297 {
<a name="l00298"></a>00298 <span class="keywordtype">double</span> sum = A[i][j];
<a name="l00299"></a>00299
<a name="l00300"></a>00300 A[j][i] = 0;
<a name="l00301"></a>00301
<a name="l00302"></a>00302 <span class="keywordflow">for</span> ( <span class="keywordtype">unsigned</span> <span class="keywordtype">int</span> k = 0; k < j; ++k )
<a name="l00303"></a>00303 sum -= A[i][k] * A[j][k];
<a name="l00304"></a>00304
<a name="l00305"></a>00305 <span class="keywordflow">if</span> (i == j)
<a name="l00306"></a>00306 {
<a name="l00307"></a>00307 <span class="keywordflow">if</span> (sum < 0) <span class="keywordflow">return</span> EXIT_FAILURE;
<a name="l00308"></a>00308
<a name="l00309"></a>00309 sum = sqrt (sum);
<a name="l00310"></a>00310
<a name="l00311"></a>00311 <span class="keywordflow">if</span> ( fabs(sum) < DBL_EPSILON ) <span class="keywordflow">return</span> EXIT_FAILURE;
<a name="l00312"></a>00312
<a name="l00313"></a>00313 A[i][j] = sum;
<a name="l00314"></a>00314 }
<a name="l00315"></a>00315 <span class="keywordflow">else</span>
<a name="l00316"></a>00316 A[i][j] = sum / A[j][j];
<a name="l00317"></a>00317 }
<a name="l00318"></a>00318
<a name="l00319"></a>00319 <span class="keywordflow">return</span> EXIT_SUCCESS;
<a name="l00320"></a>00320 }
<a name="l00321"></a>00321
<a name="l00322"></a>00322
<a name="l00323"></a>00323 <span class="keywordtype">int</span>
<a name="l00324"></a>00324 <a class="code" href="namespacehippodraw_1_1Numeric.html#a6e34a1848b8e6e40e354f0c93e6b205d" title="Solves the equation LL&#39;x = b where L is lower triangular matrix.">cholBackSolve</a> ( <span class="keyword">const</span> std::vector < std::vector< double > > & L,
<a name="l00325"></a>00325 std::vector< double > & x,
<a name="l00326"></a>00326 <span class="keyword">const</span> std::vector< double > & b )
<a name="l00327"></a>00327 {
<a name="l00328"></a>00328 <span class="keywordtype">unsigned</span> <span class="keywordtype">int</span> n = L.size();
<a name="l00329"></a>00329
<a name="l00330"></a>00330 <span class="keywordtype">double</span> sum;
<a name="l00331"></a>00331
<a name="l00332"></a>00332 <span class="comment">// Solving L y = b </span>
<a name="l00333"></a>00333 <span class="keywordflow">for</span> ( <span class="keywordtype">unsigned</span> <span class="keywordtype">int</span> i = 0; i < n ; i++)
<a name="l00334"></a>00334 {
<a name="l00335"></a>00335 sum = b [i];
<a name="l00336"></a>00336 <span class="keywordflow">for</span> ( <span class="keywordtype">unsigned</span> <span class="keywordtype">int</span> j = 0; j < i ; ++j)
<a name="l00337"></a>00337 sum -= x[j] * L[i][j];
<a name="l00338"></a>00338 x[i] = sum / L[i][i];
<a name="l00339"></a>00339 }
<a name="l00340"></a>00340
<a name="l00341"></a>00341 <span class="comment">// Solving L' x = y</span>
<a name="l00342"></a>00342 <span class="keywordflow">for</span> ( <span class="keywordtype">int</span> i = n - 1; i >= 0; i-- )
<a name="l00343"></a>00343 {
<a name="l00344"></a>00344 sum = x[i];
<a name="l00345"></a>00345 <span class="keywordflow">for</span> ( <span class="keywordtype">unsigned</span> <span class="keywordtype">int</span> j = i + 1; j < n ; j++)
<a name="l00346"></a>00346 sum -= x[j] * L[j][i];
<a name="l00347"></a>00347 x[i] = sum / L[i][i];
<a name="l00348"></a>00348 }
<a name="l00349"></a>00349
<a name="l00350"></a>00350 <span class="keywordflow">return</span> EXIT_SUCCESS;
<a name="l00351"></a>00351 }
<a name="l00352"></a>00352
<a name="l00353"></a>00353
<a name="l00354"></a>00354 <span class="keywordtype">int</span>
<a name="l00355"></a>00355 <a class="code" href="namespacehippodraw_1_1Numeric.html#a87403e87bbdd901c7e42f4bf6f314d5d" title="Inverts a SPD matrix a to get inverse Ainv using the cholesky factorization.">invertMatrix</a> ( <span class="keyword">const</span> std::vector < std::vector< double > > & A,
<a name="l00356"></a>00356 std::vector < std::vector < double > > & Ainv )
<a name="l00357"></a>00357 {
<a name="l00358"></a>00358 <span class="keywordtype">unsigned</span> <span class="keywordtype">int</span> n = A.size();
<a name="l00359"></a>00359 vector< double > x, b;
<a name="l00360"></a>00360 vector< vector< double > > L;
<a name="l00361"></a>00361
<a name="l00362"></a>00362 <span class="comment">// Set appropriate sizes for the vectors and matrices </span>
<a name="l00363"></a>00363 x.resize( n, 0.0 );
<a name="l00364"></a>00364 b.resize( n, 0.0 );
<a name="l00365"></a>00365
<a name="l00366"></a>00366 L.resize( n );
<a name="l00367"></a>00367 Ainv.resize( n );
<a name="l00368"></a>00368
<a name="l00369"></a>00369 <span class="keywordflow">for</span> ( <span class="keywordtype">unsigned</span> <span class="keywordtype">int</span> i = 0; i < n; i++ )
<a name="l00370"></a>00370 {
<a name="l00371"></a>00371 L[i].clear ();
<a name="l00372"></a>00372 L[i].resize ( n, 0.0 );
<a name="l00373"></a>00373
<a name="l00374"></a>00374 Ainv[i].clear ();
<a name="l00375"></a>00375 Ainv[i].resize ( n, 0.0 );
<a name="l00376"></a>00376 }
<a name="l00377"></a>00377
<a name="l00378"></a>00378 <span class="keywordflow">for</span> ( <span class="keywordtype">unsigned</span> <span class="keywordtype">int</span> i = 0; i < n; i++ )
<a name="l00379"></a>00379 <span class="keywordflow">for</span> ( <span class="keywordtype">unsigned</span> <span class="keywordtype">int</span> j = 0; j < n; j++ )
<a name="l00380"></a>00380 L[i][j] = A[i][j];
<a name="l00381"></a>00381
<a name="l00382"></a>00382 <span class="comment">// Do a cholesky factorization writing the lower triangular factor</span>
<a name="l00383"></a>00383 <span class="comment">// in the lower triangular part of L and setting the upper part as 0</span>
<a name="l00384"></a>00384 <a class="code" href="namespacehippodraw_1_1Numeric.html#aecbdf19ba7109e37b7e89e490aaa453a" title="The subroutine which does cholesky factorization of a given Symmetric positive definite...">cholFactor</a>( L );
<a name="l00385"></a>00385
<a name="l00386"></a>00386 <span class="keywordflow">for</span> ( <span class="keywordtype">unsigned</span> <span class="keywordtype">int</span> j = 0; j < n; j++ )
<a name="l00387"></a>00387 {
<a name="l00388"></a>00388 <span class="comment">// Set up right hand side i.e. set b = ei </span>
<a name="l00389"></a>00389 <span class="keywordflow">for</span> ( <span class="keywordtype">unsigned</span> <span class="keywordtype">int</span> i = 0; i < n; i++ )
<a name="l00390"></a>00390 b[i] = ( i == j ) ? 1 : 0;
<a name="l00391"></a>00391
<a name="l00392"></a>00392 <span class="comment">// LL'x= b is being solved here</span>
<a name="l00393"></a>00393 <a class="code" href="namespacehippodraw_1_1Numeric.html#a6e34a1848b8e6e40e354f0c93e6b205d" title="Solves the equation LL&#39;x = b where L is lower triangular matrix.">cholBackSolve</a>( L, x, b );
<a name="l00394"></a>00394
<a name="l00395"></a>00395 <span class="comment">// This x constitutes the j th coloumn of the inverse</span>
<a name="l00396"></a>00396 <span class="keywordflow">for</span> ( <span class="keywordtype">unsigned</span> <span class="keywordtype">int</span> i = 0; i < n; i++ )
<a name="l00397"></a>00397 Ainv[i][j] = x[i] ;
<a name="l00398"></a>00398 }
<a name="l00399"></a>00399
<a name="l00400"></a>00400 <span class="keywordflow">return</span> EXIT_SUCCESS;
<a name="l00401"></a>00401 }
<a name="l00402"></a>00402
<a name="l00403"></a>00403 <span class="keywordtype">int</span>
<a name="l00404"></a>00404 <a class="code" href="namespacehippodraw_1_1Numeric.html#a3c1c282c14b12b9f982f4f3f359fde67" title="Creates an n x n identity matrix for M.">eye</a> ( std::vector < std::vector < double > >& I, <span class="keywordtype">int</span> n )
<a name="l00405"></a>00405 {
<a name="l00406"></a>00406 I.resize( n );
<a name="l00407"></a>00407 <span class="keywordflow">for</span>( <span class="keywordtype">int</span> i = 0; i < n; i++ )
<a name="l00408"></a>00408 {
<a name="l00409"></a>00409 I[i].clear();
<a name="l00410"></a>00410 I[i].resize( n, 0.0 );
<a name="l00411"></a>00411 }
<a name="l00412"></a>00412
<a name="l00413"></a>00413 <span class="keywordflow">for</span>( <span class="keywordtype">int</span> i = 0; i < n; i++ )
<a name="l00414"></a>00414 I[i][i] = 1.0;
<a name="l00415"></a>00415
<a name="l00416"></a>00416 <span class="keywordflow">return</span> EXIT_SUCCESS;
<a name="l00417"></a>00417 }
<a name="l00418"></a>00418
<a name="l00419"></a>00419 <span class="keywordtype">int</span>
<a name="l00420"></a>00420 <a class="code" href="namespacehippodraw_1_1Numeric.html#a339393d88f1319ad9ae63afad77ec304" title="Given the vector it writes it to std stream.">write</a> ( <span class="keyword">const</span> std::vector < double > & a )
<a name="l00421"></a>00421 {
<a name="l00422"></a>00422 <span class="keywordtype">unsigned</span> <span class="keywordtype">int</span> n = a.size();
<a name="l00423"></a>00423
<a name="l00424"></a>00424 <span class="keywordflow">for</span> ( <span class="keywordtype">unsigned</span> <span class="keywordtype">int</span> i = 0; i < n; ++i )
<a name="l00425"></a>00425 cout << setprecision( 15 ) << a[i] << endl;
<a name="l00426"></a>00426 cout << endl;
<a name="l00427"></a>00427
<a name="l00428"></a>00428 <span class="keywordflow">return</span> EXIT_SUCCESS;
<a name="l00429"></a>00429 }
<a name="l00430"></a>00430
<a name="l00431"></a>00431
<a name="l00432"></a>00432 <span class="keywordtype">int</span>
<a name="l00433"></a>00433 <a class="code" href="namespacehippodraw_1_1Numeric.html#a339393d88f1319ad9ae63afad77ec304" title="Given the vector it writes it to std stream.">write</a> ( <span class="keyword">const</span> std::vector < std::vector < double > > & A )
<a name="l00434"></a>00434 {
<a name="l00435"></a>00435 <span class="keywordtype">unsigned</span> <span class="keywordtype">int</span> n = A.size();
<a name="l00436"></a>00436
<a name="l00437"></a>00437 <span class="keywordflow">for</span> ( <span class="keywordtype">unsigned</span> <span class="keywordtype">int</span> i = 0; i < n; ++i )
<a name="l00438"></a>00438 {
<a name="l00439"></a>00439 <span class="keywordflow">for</span> ( <span class="keywordtype">unsigned</span> <span class="keywordtype">int</span> j = 0; j < n; ++j )
<a name="l00440"></a>00440 cout << setw( 8 ) << setprecision( 4 ) << A[i][j];
<a name="l00441"></a>00441 cout << endl;
<a name="l00442"></a>00442 }
<a name="l00443"></a>00443
<a name="l00444"></a>00444 cout << endl;
<a name="l00445"></a>00445
<a name="l00446"></a>00446 <span class="keywordflow">return</span> EXIT_SUCCESS;
<a name="l00447"></a>00447 }
<a name="l00448"></a>00448
<a name="l00449"></a>00449
<a name="l00450"></a>00450 <span class="keywordtype">int</span>
<a name="l00451"></a>00451 <a class="code" href="namespacehippodraw_1_1Numeric.html#af167b6f41b28d68a0af3813267f70be5" title="Allocates a matrix of size m x n.">allocateMatrix</a> ( std::vector < std::vector < double > > & A,
<a name="l00452"></a>00452 <span class="keywordtype">int</span> m, <span class="keywordtype">int</span> n )
<a name="l00453"></a>00453 {
<a name="l00454"></a>00454 A.resize( m );
<a name="l00455"></a>00455 <span class="keywordflow">for</span>( <span class="keywordtype">int</span> i = 0; i < m; i++ )
<a name="l00456"></a>00456 {
<a name="l00457"></a>00457 A[i].clear();
<a name="l00458"></a>00458 A[i].resize( n, 0.0 );
<a name="l00459"></a>00459 }
<a name="l00460"></a>00460
<a name="l00461"></a>00461 <span class="keywordflow">return</span> EXIT_SUCCESS;
<a name="l00462"></a>00462 }
<a name="l00463"></a>00463
<a name="l00464"></a>00464
<a name="l00465"></a>00465 <span class="keywordtype">int</span>
<a name="l00466"></a>00466 <a class="code" href="namespacehippodraw_1_1Numeric.html#a9958ce4dea9d284442d0e5239793df76" title="Allocates a vector of size n.">allocateVector</a> ( std::vector < double > & x, <span class="keywordtype">int</span> n )
<a name="l00467"></a>00467 {
<a name="l00468"></a>00468 x.clear();
<a name="l00469"></a>00469 x.resize( n, 0.0 );
<a name="l00470"></a>00470
<a name="l00471"></a>00471 <span class="keywordflow">return</span> EXIT_SUCCESS;
<a name="l00472"></a>00472 }
<a name="l00473"></a>00473
<a name="l00474"></a>00474
<a name="l00475"></a>00475 <span class="comment">/* The driver main subroutine to check a few of above function.</span>
<a name="l00476"></a>00476 <span class="comment"> Test with the following matlab code:</span>
<a name="l00477"></a>00477 <span class="comment"> n = 4; L = tril( randn(n,n), -1 ) + eye(n,n); A = L' * L; B = randn(n , n); x = randn(n , 1); y = randn(n , 1); save test.txt -ascii A B x y;</span>
<a name="l00478"></a>00478 <span class="comment"> Then perform the operations as stated in the following program.</span>
<a name="l00479"></a>00479 <span class="comment">*/</span>
<a name="l00480"></a>00480 <span class="comment">/*int main()</span>
<a name="l00481"></a>00481 <span class="comment">{</span>
<a name="l00482"></a>00482 <span class="comment"> int n = 4;</span>
<a name="l00483"></a>00483 <span class="comment"> vector< double > x, y, z;</span>
<a name="l00484"></a>00484 <span class="comment"> vector< vector < double > > A, B, Ainv;</span>
<a name="l00485"></a>00485 <span class="comment"> ifstream fin( "test.txt" );</span>
<a name="l00486"></a>00486 <span class="comment"></span>
<a name="l00487"></a>00487 <span class="comment"> if( !fin )</span>
<a name="l00488"></a>00488 <span class="comment"> {</span>
<a name="l00489"></a>00489 <span class="comment"> cout << " Could not open the file for reading " << endl;</span>
<a name="l00490"></a>00490 <span class="comment"> return EXIT_SUCCESS;</span>
<a name="l00491"></a>00491 <span class="comment"> }</span>
<a name="l00492"></a>00492 <span class="comment"></span>
<a name="l00493"></a>00493 <span class="comment"> allocateMatrix( A, n, n );</span>
<a name="l00494"></a>00494 <span class="comment"> allocateMatrix( B, n, n );</span>
<a name="l00495"></a>00495 <span class="comment"> allocateVector( x, n );</span>
<a name="l00496"></a>00496 <span class="comment"> allocateVector( y, n );</span>
<a name="l00497"></a>00497 <span class="comment"> </span>
<a name="l00498"></a>00498 <span class="comment"> for( int i = 0; i < n; i++ )</span>
<a name="l00499"></a>00499 <span class="comment"> for( int j = 0; j < n; j++ )</span>
<a name="l00500"></a>00500 <span class="comment"> fin >> A[i][j];</span>
<a name="l00501"></a>00501 <span class="comment"> </span>
<a name="l00502"></a>00502 <span class="comment"> for( int i = 0; i < n; i++ )</span>
<a name="l00503"></a>00503 <span class="comment"> for( int j = 0; j < n; j++ )</span>
<a name="l00504"></a>00504 <span class="comment"> fin >> B[i][j];</span>
<a name="l00505"></a>00505 <span class="comment"> </span>
<a name="l00506"></a>00506 <span class="comment"> for( int i = 0; i < n; i++ )</span>
<a name="l00507"></a>00507 <span class="comment"> fin >> x[i];</span>
<a name="l00508"></a>00508 <span class="comment"> </span>
<a name="l00509"></a>00509 <span class="comment"> for( int i = 0; i < n; i++ )</span>
<a name="l00510"></a>00510 <span class="comment"> fin >> y[i];</span>
<a name="l00511"></a>00511 <span class="comment"> </span>
<a name="l00512"></a>00512 <span class="comment"> fin.close();</span>
<a name="l00513"></a>00513 <span class="comment"> </span>
<a name="l00514"></a>00514 <span class="comment"> //cout << "x + 3.14 * x + y /1.4141 - 2.8 * A * x - y + (y' * A)'= " << endl;</span>
<a name="l00515"></a>00515 <span class="comment"> //write( x + 3.141 * x + y / 1.4141 - 2.8 * A * x - y + y * A );</span>
<a name="l00516"></a>00516 <span class="comment"> </span>
<a name="l00517"></a>00517 <span class="comment"> //cout << "A + 3.14 * B + A /1.4141 - A * B + 65.0 * x * y' - B = " << endl;</span>
<a name="l00518"></a>00518 <span class="comment"> //write( A + 3.14 * B + A /1.4141 - A * B + 65.0 * x * y ) - B );</span>
<a name="l00519"></a>00519 <span class="comment"></span>
<a name="l00520"></a>00520 <span class="comment"> //cout << "inv(A) = " << endl;</span>
<a name="l00521"></a>00521 <span class="comment"> //invertMatrix( A, Ainv );</span>
<a name="l00522"></a>00522 <span class="comment"> //write( Ainv );</span>
<a name="l00523"></a>00523 <span class="comment"></span>
<a name="l00524"></a>00524 <span class="comment"> </span>
<a name="l00525"></a>00525 <span class="comment"> double temp = ( 1 + innerProduct( y, B * y ) / ys ) / ys;</span>
<a name="l00526"></a>00526 <span class="comment"> </span>
<a name="l00527"></a>00527 <span class="comment"> t1 = ( outerProduct(s, y) * B)/ys + ( m_M * outerProduct(y , s) ) / ys;</span>
<a name="l00528"></a>00528 <span class="comment"> t2 = temp * outerProduct( s, s );</span>
<a name="l00529"></a>00529 <span class="comment"> B = B - t1 + t2;</span>
<a name="l00530"></a>00530 <span class="comment"> </span>
<a name="l00531"></a>00531 <span class="comment"> </span>
<a name="l00532"></a>00532 <span class="comment"> return EXIT_SUCCESS;</span>
<a name="l00533"></a>00533 <span class="comment">}</span>
<a name="l00534"></a>00534 <span class="comment">*/</span>
<a name="l00535"></a>00535
<a name="l00536"></a>00536 } <span class="comment">// namespace Numeric</span>
<a name="l00537"></a>00537 } <span class="comment">// namespace hippodraw</span>
</pre></div></div>
<hr size="1"><address style="align: right;"><small>
Generated for HippoDraw Class Library by <a href="http://www.doxygen.org/index.html">
<img src="doxygen.png" alt="doxygen" align="middle" border=0>
</a> </small></address>
</body>
</html>
