<!DOCTYPE html PUBLIC "-//W3C//DTD XHTML 1.0 Transitional//EN" "http://www.w3.org/TR/xhtml1/DTD/xhtml1-transitional.dtd">
<html xmlns="http://www.w3.org/1999/xhtml">
<head>
<meta http-equiv="Content-Type" content="text/xhtml;charset=UTF-8"/>
<meta http-equiv="X-UA-Compatible" content="IE=9"/>
<title>CPPLapack: dgematrix-lapack.hpp Source File</title>

<link href="tabs.css" rel="stylesheet" type="text/css"/>
<link href="doxygen.css" rel="stylesheet" type="text/css" />

<link href="search/search.css" rel="stylesheet" type="text/css"/>
<script type="text/javascript" src="jquery.js"></script>
<script type="text/javascript" src="search/search.js"></script>
<script type="text/javascript">
  $(document).ready(function() { searchBox.OnSelectItem(0); });
</script>

</head>
<body>
<div id="top"><!-- do not remove this div! -->


<div id="titlearea">
<table cellspacing="0" cellpadding="0">
 <tbody>
 <tr style="height: 56px;">
  
  
  <td style="padding-left: 0.5em;">
   <div id="projectname">CPPLapack
   
   </div>
   
  </td>
  
  
  
 </tr>
 </tbody>
</table>
</div>

<!-- Generated by Doxygen 1.7.6.1 -->
<script type="text/javascript">
var searchBox = new SearchBox("searchBox", "search",false,'Search');
</script>
  <div id="navrow1" class="tabs">
    <ul class="tablist">
      <li><a href="index.html"><span>Main&#160;Page</span></a></li>
      <li><a href="pages.html"><span>Related&#160;Pages</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>
      <li>
        <div id="MSearchBox" class="MSearchBoxInactive">
        <span class="left">
          <img id="MSearchSelect" src="search/mag_sel.png"
               onmouseover="return searchBox.OnSearchSelectShow()"
               onmouseout="return searchBox.OnSearchSelectHide()"
               alt=""/>
          <input type="text" id="MSearchField" value="Search" accesskey="S"
               onfocus="searchBox.OnSearchFieldFocus(true)" 
               onblur="searchBox.OnSearchFieldFocus(false)" 
               onkeyup="searchBox.OnSearchFieldChange(event)"/>
          </span><span class="right">
            <a id="MSearchClose" href="javascript:searchBox.CloseResultsWindow()"><img id="MSearchCloseImg" border="0" src="search/close.png" alt=""/></a>
          </span>
        </div>
      </li>
    </ul>
  </div>
  <div id="navrow2" class="tabs2">
    <ul class="tablist">
      <li><a href="files.html"><span>File&#160;List</span></a></li>
      <li><a href="globals.html"><span>File&#160;Members</span></a></li>
    </ul>
  </div>
<!-- window showing the filter options -->
<div id="MSearchSelectWindow"
     onmouseover="return searchBox.OnSearchSelectShow()"
     onmouseout="return searchBox.OnSearchSelectHide()"
     onkeydown="return searchBox.OnSearchSelectKey(event)">
<a class="SelectItem" href="javascript:void(0)" onclick="searchBox.OnSelectItem(0)"><span class="SelectionMark">&#160;</span>All</a><a class="SelectItem" href="javascript:void(0)" onclick="searchBox.OnSelectItem(1)"><span class="SelectionMark">&#160;</span>Classes</a><a class="SelectItem" href="javascript:void(0)" onclick="searchBox.OnSelectItem(2)"><span class="SelectionMark">&#160;</span>Files</a><a class="SelectItem" href="javascript:void(0)" onclick="searchBox.OnSelectItem(3)"><span class="SelectionMark">&#160;</span>Functions</a><a class="SelectItem" href="javascript:void(0)" onclick="searchBox.OnSelectItem(4)"><span class="SelectionMark">&#160;</span>Variables</a><a class="SelectItem" href="javascript:void(0)" onclick="searchBox.OnSelectItem(5)"><span class="SelectionMark">&#160;</span>Friends</a></div>

<!-- iframe showing the search results (closed by default) -->
<div id="MSearchResultsWindow">
<iframe src="javascript:void(0)" frameborder="0" 
        name="MSearchResults" id="MSearchResults">
</iframe>
</div>

  <div id="nav-path" class="navpath">
    <ul>
      <li class="navelem"><a class="el" href="dir_a37d0dd275f7ec49ee26f1d33d1059f5.html">include</a>      </li>
      <li class="navelem"><a class="el" href="dir_fd252eee318a6d389ebe24aa4cc58da7.html">dgematrix-</a>      </li>
    </ul>
  </div>
</div>
<div class="header">
  <div class="headertitle">
<div class="title">dgematrix-lapack.hpp</div>  </div>
</div><!--header-->
<div class="contents">
<a href="dgematrix-lapack_8hpp.html">Go to the documentation of this file.</a><div class="fragment"><pre class="fragment"><a name="l00001"></a>00001 <span class="comment">//=============================================================================</span><span class="comment"></span>
<a name="l00002"></a>00002 <span class="comment">/*! solve A*X=Y using dgesv\n</span>
<a name="l00003"></a>00003 <span class="comment">  The argument is dgematrix Y. Y is overwritten and become the solution X.</span>
<a name="l00004"></a>00004 <span class="comment">  A is also overwritten and become P*l*U.</span>
<a name="l00005"></a>00005 <span class="comment">*/</span>
<a name="l00006"></a><a class="code" href="classdgematrix.html#ae25bb3ac4f085be249786264c3a0328e">00006</a> <span class="keyword">inline</span> <span class="keywordtype">long</span> <a class="code" href="classdgematrix.html#ae25bb3ac4f085be249786264c3a0328e">dgematrix::dgesv</a>(<a class="code" href="classdgematrix.html" title="Real Double-precision General Dence Matrix Class.">dgematrix</a>&amp; mat)
<a name="l00007"></a>00007 {VERBOSE_REPORT;
<a name="l00008"></a>00008 <span class="preprocessor">#ifdef  CPPL_DEBUG</span>
<a name="l00009"></a>00009 <span class="preprocessor"></span>  <span class="keywordflow">if</span>(<a class="code" href="classdgematrix.html#a3852aa1d0b876a74c1070e982cf388a2" title="matrix row size">m</a>!=<a class="code" href="classdgematrix.html#a60d26e3a499a235310e59a0aebbe0d0a" title="matrix column size">n</a> || <a class="code" href="classdgematrix.html#a3852aa1d0b876a74c1070e982cf388a2" title="matrix row size">m</a>!=mat.<a class="code" href="classdgematrix.html#a3852aa1d0b876a74c1070e982cf388a2" title="matrix row size">m</a>){
<a name="l00010"></a>00010     ERROR_REPORT;
<a name="l00011"></a>00011     std::cerr &lt;&lt; <span class="stringliteral">&quot;These two matrices cannot be solved.&quot;</span> &lt;&lt; std::endl
<a name="l00012"></a>00012               &lt;&lt; <span class="stringliteral">&quot;Your input was (&quot;</span> &lt;&lt; <a class="code" href="classdgematrix.html#a3852aa1d0b876a74c1070e982cf388a2" title="matrix row size">m</a> &lt;&lt; <span class="stringliteral">&quot;x&quot;</span> &lt;&lt; <a class="code" href="classdgematrix.html#a60d26e3a499a235310e59a0aebbe0d0a" title="matrix column size">n</a> &lt;&lt; <span class="stringliteral">&quot;) and (&quot;</span> &lt;&lt; mat.<a class="code" href="classdgematrix.html#a3852aa1d0b876a74c1070e982cf388a2" title="matrix row size">m</a> &lt;&lt; <span class="stringliteral">&quot;x&quot;</span> &lt;&lt; mat.<a class="code" href="classdgematrix.html#a60d26e3a499a235310e59a0aebbe0d0a" title="matrix column size">n</a> &lt;&lt; <span class="stringliteral">&quot;).&quot;</span> &lt;&lt; std::endl;
<a name="l00013"></a>00013     exit(1);
<a name="l00014"></a>00014   }
<a name="l00015"></a>00015 <span class="preprocessor">#endif//CPPL_DEBUG </span>
<a name="l00016"></a>00016 <span class="preprocessor"></span>  
<a name="l00017"></a>00017   <span class="keywordtype">long</span> NRHS(mat.<a class="code" href="classdgematrix.html#a60d26e3a499a235310e59a0aebbe0d0a" title="matrix column size">n</a>), LDA(<a class="code" href="classdgematrix.html#a60d26e3a499a235310e59a0aebbe0d0a" title="matrix column size">n</a>), *IPIV(<span class="keyword">new</span> <span class="keywordtype">long</span>[<a class="code" href="classdgematrix.html#a60d26e3a499a235310e59a0aebbe0d0a" title="matrix column size">n</a>]), LDB(mat.<a class="code" href="classdgematrix.html#a3852aa1d0b876a74c1070e982cf388a2" title="matrix row size">m</a>), INFO(1);
<a name="l00018"></a>00018   dgesv_(n, NRHS, <a class="code" href="classdgematrix.html#a006ce8c2f9add380bcb2816fdb13aa01" title="1D array to store matrix data">array</a>, LDA, IPIV, mat.<a class="code" href="classdgematrix.html#a006ce8c2f9add380bcb2816fdb13aa01" title="1D array to store matrix data">array</a>, LDB, INFO);
<a name="l00019"></a>00019   <span class="keyword">delete</span> [] IPIV;
<a name="l00020"></a>00020   
<a name="l00021"></a>00021   <span class="keywordflow">if</span>(INFO!=0){
<a name="l00022"></a>00022     WARNING_REPORT;
<a name="l00023"></a>00023     std::cerr &lt;&lt; <span class="stringliteral">&quot;Serious trouble happend. INFO = &quot;</span> &lt;&lt; INFO &lt;&lt; <span class="stringliteral">&quot;.&quot;</span> &lt;&lt; std::endl;
<a name="l00024"></a>00024   }
<a name="l00025"></a>00025   <span class="keywordflow">return</span> INFO;
<a name="l00026"></a>00026 }
<a name="l00027"></a>00027 
<a name="l00028"></a>00028 <span class="comment">//=============================================================================</span><span class="comment"></span>
<a name="l00029"></a>00029 <span class="comment">/*! solve A*x=y using dgesv\n</span>
<a name="l00030"></a>00030 <span class="comment">  The argument is dcovector y. y is overwritten and become the solution x.</span>
<a name="l00031"></a>00031 <span class="comment">  A is also overwritten and become P*l*U.</span>
<a name="l00032"></a>00032 <span class="comment">*/</span>
<a name="l00033"></a><a class="code" href="classdgematrix.html#a0ff6861efc9d6223c61184cd0026d8e9">00033</a> <span class="keyword">inline</span> <span class="keywordtype">long</span> <a class="code" href="classdgematrix.html#ae25bb3ac4f085be249786264c3a0328e">dgematrix::dgesv</a>(<a class="code" href="classdcovector.html" title="Real Double-precision Column Vector Class.">dcovector</a>&amp; vec)
<a name="l00034"></a>00034 {VERBOSE_REPORT;
<a name="l00035"></a>00035 <span class="preprocessor">#ifdef  CPPL_DEBUG</span>
<a name="l00036"></a>00036 <span class="preprocessor"></span>  <span class="keywordflow">if</span>(<a class="code" href="classdgematrix.html#a3852aa1d0b876a74c1070e982cf388a2" title="matrix row size">m</a>!=<a class="code" href="classdgematrix.html#a60d26e3a499a235310e59a0aebbe0d0a" title="matrix column size">n</a> || <a class="code" href="classdgematrix.html#a3852aa1d0b876a74c1070e982cf388a2" title="matrix row size">m</a>!=vec.<a class="code" href="classdcovector.html#a08e64e5ecb5121100787d0f5d48f0af4" title="vector size">l</a>){
<a name="l00037"></a>00037     ERROR_REPORT;
<a name="l00038"></a>00038     std::cerr &lt;&lt; <span class="stringliteral">&quot;These matrix and vector cannot be solved.&quot;</span> &lt;&lt; std::endl
<a name="l00039"></a>00039               &lt;&lt; <span class="stringliteral">&quot;Your input was (&quot;</span> &lt;&lt; <a class="code" href="classdgematrix.html#a3852aa1d0b876a74c1070e982cf388a2" title="matrix row size">m</a> &lt;&lt; <span class="stringliteral">&quot;x&quot;</span> &lt;&lt; <a class="code" href="classdgematrix.html#a60d26e3a499a235310e59a0aebbe0d0a" title="matrix column size">n</a> &lt;&lt; <span class="stringliteral">&quot;) and (&quot;</span> &lt;&lt; vec.<a class="code" href="classdcovector.html#a08e64e5ecb5121100787d0f5d48f0af4" title="vector size">l</a> &lt;&lt; <span class="stringliteral">&quot;).&quot;</span> &lt;&lt; std::endl;
<a name="l00040"></a>00040     exit(1);
<a name="l00041"></a>00041   }
<a name="l00042"></a>00042 <span class="preprocessor">#endif//CPPL_DEBUG </span>
<a name="l00043"></a>00043 <span class="preprocessor"></span>  
<a name="l00044"></a>00044   <span class="keywordtype">long</span> NRHS(1), LDA(<a class="code" href="classdgematrix.html#a60d26e3a499a235310e59a0aebbe0d0a" title="matrix column size">n</a>), *IPIV(<span class="keyword">new</span> <span class="keywordtype">long</span>[<a class="code" href="classdgematrix.html#a60d26e3a499a235310e59a0aebbe0d0a" title="matrix column size">n</a>]), LDB(vec.<a class="code" href="classdcovector.html#a08e64e5ecb5121100787d0f5d48f0af4" title="vector size">l</a>), INFO(1);
<a name="l00045"></a>00045   dgesv_(n, NRHS, <a class="code" href="classdgematrix.html#a006ce8c2f9add380bcb2816fdb13aa01" title="1D array to store matrix data">array</a>, LDA, IPIV, vec.<a class="code" href="classdcovector.html#a9b188bf160d8f7e24ea1d41deaceaf6d" title="1D array to store vector data">array</a>, LDB, INFO);
<a name="l00046"></a>00046   <span class="keyword">delete</span> [] IPIV;
<a name="l00047"></a>00047   
<a name="l00048"></a>00048   <span class="keywordflow">if</span>(INFO!=0){
<a name="l00049"></a>00049     WARNING_REPORT;
<a name="l00050"></a>00050     std::cerr &lt;&lt; <span class="stringliteral">&quot;Serious trouble happend. INFO = &quot;</span> &lt;&lt; INFO &lt;&lt; <span class="stringliteral">&quot;.&quot;</span> &lt;&lt; std::endl;
<a name="l00051"></a>00051   }
<a name="l00052"></a>00052   <span class="keywordflow">return</span> INFO;
<a name="l00053"></a>00053 }
<a name="l00054"></a>00054 <span class="comment"></span>
<a name="l00055"></a>00055 <span class="comment">///////////////////////////////////////////////////////////////////////////////</span>
<a name="l00056"></a>00056 <span class="comment">///////////////////////////////////////////////////////////////////////////////</span>
<a name="l00057"></a>00057 <span class="comment">///////////////////////////////////////////////////////////////////////////////</span>
<a name="l00058"></a>00058 <span class="comment"></span>
<a name="l00059"></a>00059 <span class="comment">//=============================================================================</span><span class="comment"></span>
<a name="l00060"></a>00060 <span class="comment">/*! solve overdetermined or underdetermined A*X=Y using dgels\n</span>
<a name="l00061"></a>00061 <span class="comment"> */</span>
<a name="l00062"></a><a class="code" href="classdgematrix.html#a5007de00e1868c04a553b2cbf35ed075">00062</a> <span class="keyword">inline</span> <span class="keywordtype">long</span> <a class="code" href="classdgematrix.html#a5007de00e1868c04a553b2cbf35ed075">dgematrix::dgels</a>(<a class="code" href="classdgematrix.html" title="Real Double-precision General Dence Matrix Class.">dgematrix</a>&amp; mat)
<a name="l00063"></a>00063 {VERBOSE_REPORT;
<a name="l00064"></a>00064 <span class="preprocessor">#ifdef  CPPL_DEBUG</span>
<a name="l00065"></a>00065 <span class="preprocessor"></span>  <span class="keywordflow">if</span>(<a class="code" href="classdgematrix.html#a3852aa1d0b876a74c1070e982cf388a2" title="matrix row size">m</a>!=mat.<a class="code" href="classdgematrix.html#a3852aa1d0b876a74c1070e982cf388a2" title="matrix row size">m</a>){
<a name="l00066"></a>00066     ERROR_REPORT;
<a name="l00067"></a>00067     std::cerr &lt;&lt; <span class="stringliteral">&quot;These two matrices cannot be solved.&quot;</span> &lt;&lt; std::endl
<a name="l00068"></a>00068               &lt;&lt; <span class="stringliteral">&quot;Your input was (&quot;</span> &lt;&lt; <a class="code" href="classdgematrix.html#a3852aa1d0b876a74c1070e982cf388a2" title="matrix row size">m</a> &lt;&lt; <span class="stringliteral">&quot;x&quot;</span> &lt;&lt; <a class="code" href="classdgematrix.html#a60d26e3a499a235310e59a0aebbe0d0a" title="matrix column size">n</a> &lt;&lt; <span class="stringliteral">&quot;) and (&quot;</span> &lt;&lt; mat.<a class="code" href="classdgematrix.html#a3852aa1d0b876a74c1070e982cf388a2" title="matrix row size">m</a> &lt;&lt; <span class="stringliteral">&quot;x&quot;</span> &lt;&lt; mat.<a class="code" href="classdgematrix.html#a60d26e3a499a235310e59a0aebbe0d0a" title="matrix column size">n</a> &lt;&lt; <span class="stringliteral">&quot;).&quot;</span> &lt;&lt; std::endl;
<a name="l00069"></a>00069     exit(1);
<a name="l00070"></a>00070   }
<a name="l00071"></a>00071 <span class="preprocessor">#endif//CPPL_DEBUG    </span>
<a name="l00072"></a>00072 <span class="preprocessor"></span>  
<a name="l00073"></a>00073   <span class="keywordflow">if</span>(<a class="code" href="classdgematrix.html#a3852aa1d0b876a74c1070e982cf388a2" title="matrix row size">m</a>&lt;<a class="code" href="classdgematrix.html#a60d26e3a499a235310e59a0aebbe0d0a" title="matrix column size">n</a>){ <span class="comment">//underdetermined</span>
<a name="l00074"></a>00074     <a class="code" href="classdgematrix.html" title="Real Double-precision General Dence Matrix Class.">dgematrix</a> tmp(<a class="code" href="classdgematrix.html#a60d26e3a499a235310e59a0aebbe0d0a" title="matrix column size">n</a>,mat.<a class="code" href="classdgematrix.html#a60d26e3a499a235310e59a0aebbe0d0a" title="matrix column size">n</a>);
<a name="l00075"></a>00075     <span class="keywordflow">for</span>(<span class="keywordtype">long</span> <a class="code" href="classdgematrix.html#a67f90160391aa6dbecc09b8eb6e23246">i</a>=0; <a class="code" href="classdgematrix.html#a67f90160391aa6dbecc09b8eb6e23246">i</a>&lt;mat.<a class="code" href="classdgematrix.html#a3852aa1d0b876a74c1070e982cf388a2" title="matrix row size">m</a>; <a class="code" href="classdgematrix.html#a67f90160391aa6dbecc09b8eb6e23246">i</a>++){ <span class="keywordflow">for</span>(<span class="keywordtype">long</span> j=0; j&lt;mat.<a class="code" href="classdgematrix.html#a60d26e3a499a235310e59a0aebbe0d0a" title="matrix column size">n</a>; j++){
<a name="l00076"></a>00076       tmp(<a class="code" href="classdgematrix.html#a67f90160391aa6dbecc09b8eb6e23246">i</a>,j) =mat(<a class="code" href="classdgematrix.html#a67f90160391aa6dbecc09b8eb6e23246">i</a>,j);
<a name="l00077"></a>00077     }}
<a name="l00078"></a>00078     mat.<a class="code" href="classdgematrix.html#aff82fa7ef089e6c565bfd4aaea4b9a27">clear</a>();
<a name="l00079"></a>00079     <a class="code" href="classdgematrix.html#a363aa32a0be4faacc1427b16aa59408c">swap</a>(mat,tmp);
<a name="l00080"></a>00080   }
<a name="l00081"></a>00081   
<a name="l00082"></a>00082   <span class="keywordtype">char</span> TRANS(<span class="charliteral">&#39;n&#39;</span>);
<a name="l00083"></a>00083   <span class="keywordtype">long</span> NRHS(mat.<a class="code" href="classdgematrix.html#a60d26e3a499a235310e59a0aebbe0d0a" title="matrix column size">n</a>), LDA(<a class="code" href="classdgematrix.html#a3852aa1d0b876a74c1070e982cf388a2" title="matrix row size">m</a>), LDB(mat.<a class="code" href="classdgematrix.html#a3852aa1d0b876a74c1070e982cf388a2" title="matrix row size">m</a>),
<a name="l00084"></a>00084     LWORK(std::min(<a class="code" href="classdgematrix.html#a3852aa1d0b876a74c1070e982cf388a2" title="matrix row size">m</a>,<a class="code" href="classdgematrix.html#a60d26e3a499a235310e59a0aebbe0d0a" title="matrix column size">n</a>)+std::max(std::min(<a class="code" href="classdgematrix.html#a3852aa1d0b876a74c1070e982cf388a2" title="matrix row size">m</a>,<a class="code" href="classdgematrix.html#a60d26e3a499a235310e59a0aebbe0d0a" title="matrix column size">n</a>),NRHS)), INFO(1);
<a name="l00085"></a>00085   <span class="keywordtype">double</span> *WORK(<span class="keyword">new</span> <span class="keywordtype">double</span>[LWORK]);
<a name="l00086"></a>00086   dgels_(TRANS, <a class="code" href="classdgematrix.html#a3852aa1d0b876a74c1070e982cf388a2" title="matrix row size">m</a>, <a class="code" href="classdgematrix.html#a60d26e3a499a235310e59a0aebbe0d0a" title="matrix column size">n</a>, NRHS, <a class="code" href="classdgematrix.html#a006ce8c2f9add380bcb2816fdb13aa01" title="1D array to store matrix data">array</a>, LDA, mat.<a class="code" href="classdgematrix.html#a006ce8c2f9add380bcb2816fdb13aa01" title="1D array to store matrix data">array</a>, LDB, WORK, LWORK, INFO);
<a name="l00087"></a>00087   <span class="keyword">delete</span> [] WORK;
<a name="l00088"></a>00088   
<a name="l00089"></a>00089   <span class="keywordflow">if</span>(<a class="code" href="classdgematrix.html#a3852aa1d0b876a74c1070e982cf388a2" title="matrix row size">m</a>&gt;<a class="code" href="classdgematrix.html#a60d26e3a499a235310e59a0aebbe0d0a" title="matrix column size">n</a>){ <span class="comment">//overdetermined</span>
<a name="l00090"></a>00090     <a class="code" href="classdgematrix.html" title="Real Double-precision General Dence Matrix Class.">dgematrix</a> tmp(<a class="code" href="classdgematrix.html#a60d26e3a499a235310e59a0aebbe0d0a" title="matrix column size">n</a>,mat.<a class="code" href="classdgematrix.html#a60d26e3a499a235310e59a0aebbe0d0a" title="matrix column size">n</a>);
<a name="l00091"></a>00091     <span class="keywordflow">for</span>(<span class="keywordtype">long</span> <a class="code" href="classdgematrix.html#a67f90160391aa6dbecc09b8eb6e23246">i</a>=0; <a class="code" href="classdgematrix.html#a67f90160391aa6dbecc09b8eb6e23246">i</a>&lt;tmp.<a class="code" href="classdgematrix.html#a3852aa1d0b876a74c1070e982cf388a2" title="matrix row size">m</a>; <a class="code" href="classdgematrix.html#a67f90160391aa6dbecc09b8eb6e23246">i</a>++){ <span class="keywordflow">for</span>(<span class="keywordtype">long</span> j=0; j&lt;tmp.<a class="code" href="classdgematrix.html#a60d26e3a499a235310e59a0aebbe0d0a" title="matrix column size">n</a>; j++){
<a name="l00092"></a>00092       tmp(<a class="code" href="classdgematrix.html#a67f90160391aa6dbecc09b8eb6e23246">i</a>,j) =mat(<a class="code" href="classdgematrix.html#a67f90160391aa6dbecc09b8eb6e23246">i</a>,j);
<a name="l00093"></a>00093     }}
<a name="l00094"></a>00094     mat.<a class="code" href="classdgematrix.html#aff82fa7ef089e6c565bfd4aaea4b9a27">clear</a>();
<a name="l00095"></a>00095     <a class="code" href="classdgematrix.html#a363aa32a0be4faacc1427b16aa59408c">swap</a>(mat,tmp);
<a name="l00096"></a>00096   }
<a name="l00097"></a>00097   
<a name="l00098"></a>00098   <span class="keywordflow">if</span>(INFO!=0){
<a name="l00099"></a>00099     WARNING_REPORT;
<a name="l00100"></a>00100     std::cerr &lt;&lt; <span class="stringliteral">&quot;Serious trouble happend. INFO = &quot;</span> &lt;&lt; INFO &lt;&lt; <span class="stringliteral">&quot;.&quot;</span> &lt;&lt; std::endl;
<a name="l00101"></a>00101   }
<a name="l00102"></a>00102   <span class="keywordflow">return</span> INFO;
<a name="l00103"></a>00103 }
<a name="l00104"></a>00104 
<a name="l00105"></a>00105 <span class="comment">//=============================================================================</span><span class="comment"></span>
<a name="l00106"></a>00106 <span class="comment">/*! solve overdetermined or underdetermined A*x=y using dgels\n</span>
<a name="l00107"></a>00107 <span class="comment"> */</span>
<a name="l00108"></a><a class="code" href="classdgematrix.html#ae91ffc25042f1f25977c048953a6ad41">00108</a> <span class="keyword">inline</span> <span class="keywordtype">long</span> <a class="code" href="classdgematrix.html#a5007de00e1868c04a553b2cbf35ed075">dgematrix::dgels</a>(<a class="code" href="classdcovector.html" title="Real Double-precision Column Vector Class.">dcovector</a>&amp; vec)
<a name="l00109"></a>00109 {VERBOSE_REPORT;
<a name="l00110"></a>00110 <span class="preprocessor">#ifdef  CPPL_DEBUG</span>
<a name="l00111"></a>00111 <span class="preprocessor"></span>  <span class="keywordflow">if</span>(<a class="code" href="classdgematrix.html#a3852aa1d0b876a74c1070e982cf388a2" title="matrix row size">m</a>!=vec.<a class="code" href="classdcovector.html#a08e64e5ecb5121100787d0f5d48f0af4" title="vector size">l</a>){
<a name="l00112"></a>00112     ERROR_REPORT;
<a name="l00113"></a>00113     std::cerr &lt;&lt; <span class="stringliteral">&quot;These matrix and vector cannot be solved.&quot;</span> &lt;&lt; std::endl
<a name="l00114"></a>00114               &lt;&lt; <span class="stringliteral">&quot;Your input was (&quot;</span> &lt;&lt; <a class="code" href="classdgematrix.html#a3852aa1d0b876a74c1070e982cf388a2" title="matrix row size">m</a> &lt;&lt; <span class="stringliteral">&quot;x&quot;</span> &lt;&lt; <a class="code" href="classdgematrix.html#a60d26e3a499a235310e59a0aebbe0d0a" title="matrix column size">n</a> &lt;&lt; <span class="stringliteral">&quot;) and (&quot;</span> &lt;&lt; vec.<a class="code" href="classdcovector.html#a08e64e5ecb5121100787d0f5d48f0af4" title="vector size">l</a> &lt;&lt; <span class="stringliteral">&quot;).&quot;</span> &lt;&lt; std::endl;
<a name="l00115"></a>00115     exit(1);
<a name="l00116"></a>00116   }
<a name="l00117"></a>00117 <span class="preprocessor">#endif//CPPL_DEBUG    </span>
<a name="l00118"></a>00118 <span class="preprocessor"></span>  
<a name="l00119"></a>00119   <span class="keywordflow">if</span>(<a class="code" href="classdgematrix.html#a3852aa1d0b876a74c1070e982cf388a2" title="matrix row size">m</a>&lt;<a class="code" href="classdgematrix.html#a60d26e3a499a235310e59a0aebbe0d0a" title="matrix column size">n</a>){ <span class="comment">//underdetermined</span>
<a name="l00120"></a>00120     <a class="code" href="classdcovector.html" title="Real Double-precision Column Vector Class.">dcovector</a> tmp(<a class="code" href="classdgematrix.html#a60d26e3a499a235310e59a0aebbe0d0a" title="matrix column size">n</a>);
<a name="l00121"></a>00121     <span class="keywordflow">for</span>(<span class="keywordtype">long</span> <a class="code" href="classdgematrix.html#a67f90160391aa6dbecc09b8eb6e23246">i</a>=0; <a class="code" href="classdgematrix.html#a67f90160391aa6dbecc09b8eb6e23246">i</a>&lt;vec.<a class="code" href="classdcovector.html#a08e64e5ecb5121100787d0f5d48f0af4" title="vector size">l</a>; <a class="code" href="classdgematrix.html#a67f90160391aa6dbecc09b8eb6e23246">i</a>++){ tmp(<a class="code" href="classdgematrix.html#a67f90160391aa6dbecc09b8eb6e23246">i</a>)=vec(<a class="code" href="classdgematrix.html#a67f90160391aa6dbecc09b8eb6e23246">i</a>); }
<a name="l00122"></a>00122     vec.<a class="code" href="classdcovector.html#aa0a1e4ad0a058308b5c839609f05e973">clear</a>();
<a name="l00123"></a>00123     <a class="code" href="classdgematrix.html#a363aa32a0be4faacc1427b16aa59408c">swap</a>(vec,tmp);
<a name="l00124"></a>00124   }
<a name="l00125"></a>00125   
<a name="l00126"></a>00126   <span class="keywordtype">char</span> TRANS(<span class="charliteral">&#39;n&#39;</span>);
<a name="l00127"></a>00127   <span class="keywordtype">long</span> NRHS(1), LDA(<a class="code" href="classdgematrix.html#a3852aa1d0b876a74c1070e982cf388a2" title="matrix row size">m</a>), LDB(vec.<a class="code" href="classdcovector.html#a08e64e5ecb5121100787d0f5d48f0af4" title="vector size">l</a>),
<a name="l00128"></a>00128     LWORK(std::min(<a class="code" href="classdgematrix.html#a3852aa1d0b876a74c1070e982cf388a2" title="matrix row size">m</a>,<a class="code" href="classdgematrix.html#a60d26e3a499a235310e59a0aebbe0d0a" title="matrix column size">n</a>)+std::max(std::min(<a class="code" href="classdgematrix.html#a3852aa1d0b876a74c1070e982cf388a2" title="matrix row size">m</a>,<a class="code" href="classdgematrix.html#a60d26e3a499a235310e59a0aebbe0d0a" title="matrix column size">n</a>),NRHS)), INFO(1);
<a name="l00129"></a>00129   <span class="keywordtype">double</span> *WORK(<span class="keyword">new</span> <span class="keywordtype">double</span>[LWORK]);
<a name="l00130"></a>00130   dgels_(TRANS, <a class="code" href="classdgematrix.html#a3852aa1d0b876a74c1070e982cf388a2" title="matrix row size">m</a>, <a class="code" href="classdgematrix.html#a60d26e3a499a235310e59a0aebbe0d0a" title="matrix column size">n</a>, NRHS, <a class="code" href="classdgematrix.html#a006ce8c2f9add380bcb2816fdb13aa01" title="1D array to store matrix data">array</a>, LDA, vec.<a class="code" href="classdcovector.html#a9b188bf160d8f7e24ea1d41deaceaf6d" title="1D array to store vector data">array</a>, LDB, WORK, LWORK, INFO);
<a name="l00131"></a>00131   <span class="keyword">delete</span> [] WORK;
<a name="l00132"></a>00132   
<a name="l00133"></a>00133   <span class="keywordflow">if</span>(<a class="code" href="classdgematrix.html#a3852aa1d0b876a74c1070e982cf388a2" title="matrix row size">m</a>&gt;<a class="code" href="classdgematrix.html#a60d26e3a499a235310e59a0aebbe0d0a" title="matrix column size">n</a>){ <span class="comment">//overdetermined</span>
<a name="l00134"></a>00134     <a class="code" href="classdcovector.html" title="Real Double-precision Column Vector Class.">dcovector</a> tmp(<a class="code" href="classdgematrix.html#a60d26e3a499a235310e59a0aebbe0d0a" title="matrix column size">n</a>);
<a name="l00135"></a>00135     <span class="keywordflow">for</span>(<span class="keywordtype">long</span> <a class="code" href="classdgematrix.html#a67f90160391aa6dbecc09b8eb6e23246">i</a>=0; <a class="code" href="classdgematrix.html#a67f90160391aa6dbecc09b8eb6e23246">i</a>&lt;tmp.<a class="code" href="classdcovector.html#a08e64e5ecb5121100787d0f5d48f0af4" title="vector size">l</a>; <a class="code" href="classdgematrix.html#a67f90160391aa6dbecc09b8eb6e23246">i</a>++){ tmp(<a class="code" href="classdgematrix.html#a67f90160391aa6dbecc09b8eb6e23246">i</a>)=vec(<a class="code" href="classdgematrix.html#a67f90160391aa6dbecc09b8eb6e23246">i</a>); }
<a name="l00136"></a>00136     vec.<a class="code" href="classdcovector.html#aa0a1e4ad0a058308b5c839609f05e973">clear</a>();
<a name="l00137"></a>00137     <a class="code" href="classdgematrix.html#a363aa32a0be4faacc1427b16aa59408c">swap</a>(vec,tmp);
<a name="l00138"></a>00138   }
<a name="l00139"></a>00139   
<a name="l00140"></a>00140   <span class="keywordflow">if</span>(INFO!=0){
<a name="l00141"></a>00141     WARNING_REPORT;
<a name="l00142"></a>00142     std::cerr &lt;&lt; <span class="stringliteral">&quot;Serious trouble happend. INFO = &quot;</span> &lt;&lt; INFO &lt;&lt; <span class="stringliteral">&quot;.&quot;</span> &lt;&lt; std::endl;
<a name="l00143"></a>00143   }
<a name="l00144"></a>00144   <span class="keywordflow">return</span> INFO;
<a name="l00145"></a>00145 }
<a name="l00146"></a>00146 
<a name="l00147"></a>00147 <span class="comment">//=============================================================================</span><span class="comment"></span>
<a name="l00148"></a>00148 <span class="comment">/*! solve overdetermined or underdetermined A*X=Y using dgels</span>
<a name="l00149"></a>00149 <span class="comment">  with the sum of residual squares output\n</span>
<a name="l00150"></a>00150 <span class="comment">  The residual is set as the columnwise sum of residual squares </span>
<a name="l00151"></a>00151 <span class="comment">  for overdetermined problems</span>
<a name="l00152"></a>00152 <span class="comment">  while it is always zero for underdetermined problems.</span>
<a name="l00153"></a>00153 <span class="comment">*/</span>
<a name="l00154"></a><a class="code" href="classdgematrix.html#ab10e284667ca1e913d96f82a9f4fbd18">00154</a> <span class="keyword">inline</span> <span class="keywordtype">long</span> <a class="code" href="classdgematrix.html#a5007de00e1868c04a553b2cbf35ed075">dgematrix::dgels</a>(<a class="code" href="classdgematrix.html" title="Real Double-precision General Dence Matrix Class.">dgematrix</a>&amp; mat, <a class="code" href="classdrovector.html" title="Real Double-precision Row Vector Class.">drovector</a>&amp; residual)
<a name="l00155"></a>00155 {VERBOSE_REPORT;
<a name="l00156"></a>00156 <span class="preprocessor">#ifdef  CPPL_DEBUG</span>
<a name="l00157"></a>00157 <span class="preprocessor"></span>  <span class="keywordflow">if</span>(<a class="code" href="classdgematrix.html#a3852aa1d0b876a74c1070e982cf388a2" title="matrix row size">m</a>!=mat.<a class="code" href="classdgematrix.html#a3852aa1d0b876a74c1070e982cf388a2" title="matrix row size">m</a>){
<a name="l00158"></a>00158     ERROR_REPORT;
<a name="l00159"></a>00159     std::cerr &lt;&lt; <span class="stringliteral">&quot;These two matrices cannot be solved.&quot;</span> &lt;&lt; std::endl
<a name="l00160"></a>00160               &lt;&lt; <span class="stringliteral">&quot;Your input was (&quot;</span> &lt;&lt; <a class="code" href="classdgematrix.html#a3852aa1d0b876a74c1070e982cf388a2" title="matrix row size">m</a> &lt;&lt; <span class="stringliteral">&quot;x&quot;</span> &lt;&lt; <a class="code" href="classdgematrix.html#a60d26e3a499a235310e59a0aebbe0d0a" title="matrix column size">n</a> &lt;&lt; <span class="stringliteral">&quot;) and (&quot;</span> &lt;&lt; mat.<a class="code" href="classdgematrix.html#a3852aa1d0b876a74c1070e982cf388a2" title="matrix row size">m</a> &lt;&lt; <span class="stringliteral">&quot;x&quot;</span> &lt;&lt; mat.<a class="code" href="classdgematrix.html#a60d26e3a499a235310e59a0aebbe0d0a" title="matrix column size">n</a> &lt;&lt; <span class="stringliteral">&quot;).&quot;</span> &lt;&lt; std::endl;
<a name="l00161"></a>00161     exit(1);
<a name="l00162"></a>00162   }
<a name="l00163"></a>00163 <span class="preprocessor">#endif//CPPL_DEBUG</span>
<a name="l00164"></a>00164 <span class="preprocessor"></span>  
<a name="l00165"></a>00165   residual.<a class="code" href="classdrovector.html#a3b2c5a9ed38d2e7d42cbafc635e45db2">resize</a>(mat.<a class="code" href="classdgematrix.html#a60d26e3a499a235310e59a0aebbe0d0a" title="matrix column size">n</a>); residual.<a class="code" href="classdrovector.html#a845494ed14746353f61c962b3a5b5793">zero</a>();
<a name="l00166"></a>00166   
<a name="l00167"></a>00167   <span class="keywordflow">if</span>(<a class="code" href="classdgematrix.html#a3852aa1d0b876a74c1070e982cf388a2" title="matrix row size">m</a>&lt;<a class="code" href="classdgematrix.html#a60d26e3a499a235310e59a0aebbe0d0a" title="matrix column size">n</a>){ <span class="comment">//underdetermined</span>
<a name="l00168"></a>00168     <a class="code" href="classdgematrix.html" title="Real Double-precision General Dence Matrix Class.">dgematrix</a> tmp(<a class="code" href="classdgematrix.html#a60d26e3a499a235310e59a0aebbe0d0a" title="matrix column size">n</a>,mat.<a class="code" href="classdgematrix.html#a60d26e3a499a235310e59a0aebbe0d0a" title="matrix column size">n</a>);
<a name="l00169"></a>00169     <span class="keywordflow">for</span>(<span class="keywordtype">long</span> <a class="code" href="classdgematrix.html#a67f90160391aa6dbecc09b8eb6e23246">i</a>=0; <a class="code" href="classdgematrix.html#a67f90160391aa6dbecc09b8eb6e23246">i</a>&lt;mat.<a class="code" href="classdgematrix.html#a3852aa1d0b876a74c1070e982cf388a2" title="matrix row size">m</a>; <a class="code" href="classdgematrix.html#a67f90160391aa6dbecc09b8eb6e23246">i</a>++){ <span class="keywordflow">for</span>(<span class="keywordtype">long</span> j=0; j&lt;mat.<a class="code" href="classdgematrix.html#a60d26e3a499a235310e59a0aebbe0d0a" title="matrix column size">n</a>; j++){
<a name="l00170"></a>00170       tmp(<a class="code" href="classdgematrix.html#a67f90160391aa6dbecc09b8eb6e23246">i</a>,j) =mat(<a class="code" href="classdgematrix.html#a67f90160391aa6dbecc09b8eb6e23246">i</a>,j);
<a name="l00171"></a>00171     }}
<a name="l00172"></a>00172     mat.<a class="code" href="classdgematrix.html#aff82fa7ef089e6c565bfd4aaea4b9a27">clear</a>();
<a name="l00173"></a>00173     <a class="code" href="classdgematrix.html#a363aa32a0be4faacc1427b16aa59408c">swap</a>(mat,tmp);
<a name="l00174"></a>00174   }
<a name="l00175"></a>00175   
<a name="l00176"></a>00176   <span class="keywordtype">char</span> TRANS(<span class="charliteral">&#39;n&#39;</span>);
<a name="l00177"></a>00177   <span class="keywordtype">long</span> NRHS(mat.<a class="code" href="classdgematrix.html#a60d26e3a499a235310e59a0aebbe0d0a" title="matrix column size">n</a>), LDA(<a class="code" href="classdgematrix.html#a3852aa1d0b876a74c1070e982cf388a2" title="matrix row size">m</a>), LDB(mat.<a class="code" href="classdgematrix.html#a3852aa1d0b876a74c1070e982cf388a2" title="matrix row size">m</a>),
<a name="l00178"></a>00178     LWORK(std::min(<a class="code" href="classdgematrix.html#a3852aa1d0b876a74c1070e982cf388a2" title="matrix row size">m</a>,<a class="code" href="classdgematrix.html#a60d26e3a499a235310e59a0aebbe0d0a" title="matrix column size">n</a>)+std::max(std::min(<a class="code" href="classdgematrix.html#a3852aa1d0b876a74c1070e982cf388a2" title="matrix row size">m</a>,<a class="code" href="classdgematrix.html#a60d26e3a499a235310e59a0aebbe0d0a" title="matrix column size">n</a>),NRHS)), INFO(1);
<a name="l00179"></a>00179   <span class="keywordtype">double</span> *WORK(<span class="keyword">new</span> <span class="keywordtype">double</span>[LWORK]);
<a name="l00180"></a>00180   dgels_(TRANS, <a class="code" href="classdgematrix.html#a3852aa1d0b876a74c1070e982cf388a2" title="matrix row size">m</a>, <a class="code" href="classdgematrix.html#a60d26e3a499a235310e59a0aebbe0d0a" title="matrix column size">n</a>, NRHS, <a class="code" href="classdgematrix.html#a006ce8c2f9add380bcb2816fdb13aa01" title="1D array to store matrix data">array</a>, LDA, mat.<a class="code" href="classdgematrix.html#a006ce8c2f9add380bcb2816fdb13aa01" title="1D array to store matrix data">array</a>, LDB, WORK, LWORK, INFO);
<a name="l00181"></a>00181   <span class="keyword">delete</span> [] WORK;
<a name="l00182"></a>00182   
<a name="l00183"></a>00183   <span class="keywordflow">if</span>(<a class="code" href="classdgematrix.html#a3852aa1d0b876a74c1070e982cf388a2" title="matrix row size">m</a>&gt;<a class="code" href="classdgematrix.html#a60d26e3a499a235310e59a0aebbe0d0a" title="matrix column size">n</a>){ <span class="comment">//overdetermined</span>
<a name="l00184"></a>00184     <span class="keywordflow">for</span>(<span class="keywordtype">long</span> <a class="code" href="classdgematrix.html#a67f90160391aa6dbecc09b8eb6e23246">i</a>=0; <a class="code" href="classdgematrix.html#a67f90160391aa6dbecc09b8eb6e23246">i</a>&lt;residual.<a class="code" href="classdrovector.html#a08df95e527c85a13ed99c4a0fd2b303e" title="vector size">l</a>; <a class="code" href="classdgematrix.html#a67f90160391aa6dbecc09b8eb6e23246">i</a>++){ <span class="keywordflow">for</span>(<span class="keywordtype">long</span> j=0; j&lt;<a class="code" href="classdgematrix.html#a3852aa1d0b876a74c1070e982cf388a2" title="matrix row size">m</a>-<a class="code" href="classdgematrix.html#a60d26e3a499a235310e59a0aebbe0d0a" title="matrix column size">n</a>; j++){
<a name="l00185"></a>00185       residual(<a class="code" href="classdgematrix.html#a67f90160391aa6dbecc09b8eb6e23246">i</a>) += std::pow(mat(n+j,<a class="code" href="classdgematrix.html#a67f90160391aa6dbecc09b8eb6e23246">i</a>), 2.0);
<a name="l00186"></a>00186     }}
<a name="l00187"></a>00187     
<a name="l00188"></a>00188     <a class="code" href="classdgematrix.html" title="Real Double-precision General Dence Matrix Class.">dgematrix</a> tmp(<a class="code" href="classdgematrix.html#a60d26e3a499a235310e59a0aebbe0d0a" title="matrix column size">n</a>,mat.<a class="code" href="classdgematrix.html#a60d26e3a499a235310e59a0aebbe0d0a" title="matrix column size">n</a>);
<a name="l00189"></a>00189     <span class="keywordflow">for</span>(<span class="keywordtype">long</span> <a class="code" href="classdgematrix.html#a67f90160391aa6dbecc09b8eb6e23246">i</a>=0; <a class="code" href="classdgematrix.html#a67f90160391aa6dbecc09b8eb6e23246">i</a>&lt;tmp.<a class="code" href="classdgematrix.html#a3852aa1d0b876a74c1070e982cf388a2" title="matrix row size">m</a>; <a class="code" href="classdgematrix.html#a67f90160391aa6dbecc09b8eb6e23246">i</a>++){ <span class="keywordflow">for</span>(<span class="keywordtype">long</span> j=0; j&lt;tmp.<a class="code" href="classdgematrix.html#a60d26e3a499a235310e59a0aebbe0d0a" title="matrix column size">n</a>; j++){
<a name="l00190"></a>00190       tmp(<a class="code" href="classdgematrix.html#a67f90160391aa6dbecc09b8eb6e23246">i</a>,j) =mat(<a class="code" href="classdgematrix.html#a67f90160391aa6dbecc09b8eb6e23246">i</a>,j);
<a name="l00191"></a>00191     }}
<a name="l00192"></a>00192     mat.<a class="code" href="classdgematrix.html#aff82fa7ef089e6c565bfd4aaea4b9a27">clear</a>();
<a name="l00193"></a>00193     <a class="code" href="classdgematrix.html#a363aa32a0be4faacc1427b16aa59408c">swap</a>(mat,tmp);
<a name="l00194"></a>00194   }
<a name="l00195"></a>00195   
<a name="l00196"></a>00196   <span class="keywordflow">if</span>(INFO!=0){
<a name="l00197"></a>00197     WARNING_REPORT;
<a name="l00198"></a>00198     std::cerr &lt;&lt; <span class="stringliteral">&quot;Serious trouble happend. INFO = &quot;</span> &lt;&lt; INFO &lt;&lt; <span class="stringliteral">&quot;.&quot;</span> &lt;&lt; std::endl;
<a name="l00199"></a>00199   }
<a name="l00200"></a>00200   <span class="keywordflow">return</span> INFO;
<a name="l00201"></a>00201 }
<a name="l00202"></a>00202 
<a name="l00203"></a>00203 <span class="comment">//=============================================================================</span><span class="comment"></span>
<a name="l00204"></a>00204 <span class="comment">/*! solve overdetermined or underdetermined A*x=y using dgels</span>
<a name="l00205"></a>00205 <span class="comment">  with the sum of residual squares output\n</span>
<a name="l00206"></a>00206 <span class="comment">  The residual is set as the sum of residual squares </span>
<a name="l00207"></a>00207 <span class="comment">  for overdetermined problems</span>
<a name="l00208"></a>00208 <span class="comment">  while it is always zero for underdetermined problems.</span>
<a name="l00209"></a>00209 <span class="comment">*/</span>
<a name="l00210"></a><a class="code" href="classdgematrix.html#a63301c542bc49999259dcef0a64433ca">00210</a> <span class="keyword">inline</span> <span class="keywordtype">long</span> <a class="code" href="classdgematrix.html#a5007de00e1868c04a553b2cbf35ed075">dgematrix::dgels</a>(<a class="code" href="classdcovector.html" title="Real Double-precision Column Vector Class.">dcovector</a>&amp; vec, <span class="keywordtype">double</span>&amp; residual)
<a name="l00211"></a>00211 {VERBOSE_REPORT;
<a name="l00212"></a>00212 <span class="preprocessor">#ifdef  CPPL_DEBUG</span>
<a name="l00213"></a>00213 <span class="preprocessor"></span>  <span class="keywordflow">if</span>(<a class="code" href="classdgematrix.html#a3852aa1d0b876a74c1070e982cf388a2" title="matrix row size">m</a>!=vec.<a class="code" href="classdcovector.html#a08e64e5ecb5121100787d0f5d48f0af4" title="vector size">l</a>){
<a name="l00214"></a>00214     ERROR_REPORT;
<a name="l00215"></a>00215     std::cerr &lt;&lt; <span class="stringliteral">&quot;These matrix and vector cannot be solved.&quot;</span> &lt;&lt; std::endl
<a name="l00216"></a>00216               &lt;&lt; <span class="stringliteral">&quot;Your input was (&quot;</span> &lt;&lt; <a class="code" href="classdgematrix.html#a3852aa1d0b876a74c1070e982cf388a2" title="matrix row size">m</a> &lt;&lt; <span class="stringliteral">&quot;x&quot;</span> &lt;&lt; <a class="code" href="classdgematrix.html#a60d26e3a499a235310e59a0aebbe0d0a" title="matrix column size">n</a> &lt;&lt; <span class="stringliteral">&quot;) and (&quot;</span> &lt;&lt; vec.<a class="code" href="classdcovector.html#a08e64e5ecb5121100787d0f5d48f0af4" title="vector size">l</a> &lt;&lt; <span class="stringliteral">&quot;).&quot;</span> &lt;&lt; std::endl;
<a name="l00217"></a>00217     exit(1);
<a name="l00218"></a>00218   }
<a name="l00219"></a>00219 <span class="preprocessor">#endif//CPPL_DEBUG    </span>
<a name="l00220"></a>00220 <span class="preprocessor"></span>  
<a name="l00221"></a>00221   residual=0.0;
<a name="l00222"></a>00222   
<a name="l00223"></a>00223   <span class="keywordflow">if</span>(<a class="code" href="classdgematrix.html#a3852aa1d0b876a74c1070e982cf388a2" title="matrix row size">m</a>&lt;<a class="code" href="classdgematrix.html#a60d26e3a499a235310e59a0aebbe0d0a" title="matrix column size">n</a>){ <span class="comment">//underdetermined</span>
<a name="l00224"></a>00224     <a class="code" href="classdcovector.html" title="Real Double-precision Column Vector Class.">dcovector</a> tmp(<a class="code" href="classdgematrix.html#a60d26e3a499a235310e59a0aebbe0d0a" title="matrix column size">n</a>);
<a name="l00225"></a>00225     <span class="keywordflow">for</span>(<span class="keywordtype">long</span> <a class="code" href="classdgematrix.html#a67f90160391aa6dbecc09b8eb6e23246">i</a>=0; <a class="code" href="classdgematrix.html#a67f90160391aa6dbecc09b8eb6e23246">i</a>&lt;vec.<a class="code" href="classdcovector.html#a08e64e5ecb5121100787d0f5d48f0af4" title="vector size">l</a>; <a class="code" href="classdgematrix.html#a67f90160391aa6dbecc09b8eb6e23246">i</a>++){ tmp(<a class="code" href="classdgematrix.html#a67f90160391aa6dbecc09b8eb6e23246">i</a>)=vec(<a class="code" href="classdgematrix.html#a67f90160391aa6dbecc09b8eb6e23246">i</a>); }
<a name="l00226"></a>00226     vec.<a class="code" href="classdcovector.html#aa0a1e4ad0a058308b5c839609f05e973">clear</a>();
<a name="l00227"></a>00227     <a class="code" href="classdgematrix.html#a363aa32a0be4faacc1427b16aa59408c">swap</a>(vec,tmp);
<a name="l00228"></a>00228   }
<a name="l00229"></a>00229   
<a name="l00230"></a>00230   <span class="keywordtype">char</span> TRANS(<span class="charliteral">&#39;n&#39;</span>);
<a name="l00231"></a>00231   <span class="keywordtype">long</span> NRHS(1), LDA(<a class="code" href="classdgematrix.html#a3852aa1d0b876a74c1070e982cf388a2" title="matrix row size">m</a>), LDB(vec.<a class="code" href="classdcovector.html#a08e64e5ecb5121100787d0f5d48f0af4" title="vector size">l</a>),
<a name="l00232"></a>00232     LWORK(std::min(<a class="code" href="classdgematrix.html#a3852aa1d0b876a74c1070e982cf388a2" title="matrix row size">m</a>,<a class="code" href="classdgematrix.html#a60d26e3a499a235310e59a0aebbe0d0a" title="matrix column size">n</a>)+std::max(std::min(<a class="code" href="classdgematrix.html#a3852aa1d0b876a74c1070e982cf388a2" title="matrix row size">m</a>,<a class="code" href="classdgematrix.html#a60d26e3a499a235310e59a0aebbe0d0a" title="matrix column size">n</a>),NRHS)), INFO(1);
<a name="l00233"></a>00233   <span class="keywordtype">double</span> *WORK(<span class="keyword">new</span> <span class="keywordtype">double</span>[LWORK]);
<a name="l00234"></a>00234   dgels_(TRANS, <a class="code" href="classdgematrix.html#a3852aa1d0b876a74c1070e982cf388a2" title="matrix row size">m</a>, <a class="code" href="classdgematrix.html#a60d26e3a499a235310e59a0aebbe0d0a" title="matrix column size">n</a>, NRHS, <a class="code" href="classdgematrix.html#a006ce8c2f9add380bcb2816fdb13aa01" title="1D array to store matrix data">array</a>, LDA, vec.<a class="code" href="classdcovector.html#a9b188bf160d8f7e24ea1d41deaceaf6d" title="1D array to store vector data">array</a>, LDB, WORK, LWORK, INFO);
<a name="l00235"></a>00235   <span class="keyword">delete</span> [] WORK;
<a name="l00236"></a>00236   
<a name="l00237"></a>00237   <span class="keywordflow">if</span>(<a class="code" href="classdgematrix.html#a3852aa1d0b876a74c1070e982cf388a2" title="matrix row size">m</a>&gt;<a class="code" href="classdgematrix.html#a60d26e3a499a235310e59a0aebbe0d0a" title="matrix column size">n</a>){ <span class="comment">//overdetermined</span>
<a name="l00238"></a>00238     <span class="keywordflow">for</span>(<span class="keywordtype">long</span> <a class="code" href="classdgematrix.html#a67f90160391aa6dbecc09b8eb6e23246">i</a>=0; <a class="code" href="classdgematrix.html#a67f90160391aa6dbecc09b8eb6e23246">i</a>&lt;<a class="code" href="classdgematrix.html#a3852aa1d0b876a74c1070e982cf388a2" title="matrix row size">m</a>-<a class="code" href="classdgematrix.html#a60d26e3a499a235310e59a0aebbe0d0a" title="matrix column size">n</a>; <a class="code" href="classdgematrix.html#a67f90160391aa6dbecc09b8eb6e23246">i</a>++){ residual+=std::pow(vec(n+<a class="code" href="classdgematrix.html#a67f90160391aa6dbecc09b8eb6e23246">i</a>),2.0); }
<a name="l00239"></a>00239     
<a name="l00240"></a>00240     <a class="code" href="classdcovector.html" title="Real Double-precision Column Vector Class.">dcovector</a> tmp(n);
<a name="l00241"></a>00241     <span class="keywordflow">for</span>(<span class="keywordtype">long</span> <a class="code" href="classdgematrix.html#a67f90160391aa6dbecc09b8eb6e23246">i</a>=0; <a class="code" href="classdgematrix.html#a67f90160391aa6dbecc09b8eb6e23246">i</a>&lt;tmp.<a class="code" href="classdcovector.html#a08e64e5ecb5121100787d0f5d48f0af4" title="vector size">l</a>; <a class="code" href="classdgematrix.html#a67f90160391aa6dbecc09b8eb6e23246">i</a>++){ tmp(<a class="code" href="classdgematrix.html#a67f90160391aa6dbecc09b8eb6e23246">i</a>)=vec(<a class="code" href="classdgematrix.html#a67f90160391aa6dbecc09b8eb6e23246">i</a>); }
<a name="l00242"></a>00242     vec.<a class="code" href="classdcovector.html#aa0a1e4ad0a058308b5c839609f05e973">clear</a>();
<a name="l00243"></a>00243     <a class="code" href="classdgematrix.html#a363aa32a0be4faacc1427b16aa59408c">swap</a>(vec,tmp);
<a name="l00244"></a>00244   }
<a name="l00245"></a>00245   
<a name="l00246"></a>00246   <span class="keywordflow">if</span>(INFO!=0){
<a name="l00247"></a>00247     WARNING_REPORT;
<a name="l00248"></a>00248     std::cerr &lt;&lt; <span class="stringliteral">&quot;Serious trouble happend. INFO = &quot;</span> &lt;&lt; INFO &lt;&lt; <span class="stringliteral">&quot;.&quot;</span> &lt;&lt; std::endl;
<a name="l00249"></a>00249   }
<a name="l00250"></a>00250   <span class="keywordflow">return</span> INFO;
<a name="l00251"></a>00251 }
<a name="l00252"></a>00252 <span class="comment"></span>
<a name="l00253"></a>00253 <span class="comment">///////////////////////////////////////////////////////////////////////////////</span>
<a name="l00254"></a>00254 <span class="comment">///////////////////////////////////////////////////////////////////////////////</span>
<a name="l00255"></a>00255 <span class="comment">///////////////////////////////////////////////////////////////////////////////</span>
<a name="l00256"></a>00256 <span class="comment"></span>
<a name="l00257"></a>00257 <span class="comment">//=============================================================================</span><span class="comment"></span>
<a name="l00258"></a>00258 <span class="comment">/*! calculate the least-squares-least-norm solution for overdetermined or </span>
<a name="l00259"></a>00259 <span class="comment">  underdetermined A*x=y using dgelss\n</span>
<a name="l00260"></a>00260 <span class="comment">*/</span>
<a name="l00261"></a><a class="code" href="classdgematrix.html#ab3abb914fecbe4b373b8bab0b4395230">00261</a> <span class="keyword">inline</span> <span class="keywordtype">long</span> <a class="code" href="classdgematrix.html#ab3abb914fecbe4b373b8bab0b4395230">dgematrix::dgelss</a>(<a class="code" href="classdcovector.html" title="Real Double-precision Column Vector Class.">dcovector</a>&amp; B, <a class="code" href="classdcovector.html" title="Real Double-precision Column Vector Class.">dcovector</a>&amp; S, <span class="keywordtype">long</span>&amp; RANK,
<a name="l00262"></a>00262                               <span class="keyword">const</span> <span class="keywordtype">double</span> RCOND =-1. )
<a name="l00263"></a>00263 {VERBOSE_REPORT;
<a name="l00264"></a>00264 <span class="preprocessor">#ifdef  CPPL_DEBUG</span>
<a name="l00265"></a>00265 <span class="preprocessor"></span>  <span class="keywordflow">if</span>(<a class="code" href="classdgematrix.html#a3852aa1d0b876a74c1070e982cf388a2" title="matrix row size">m</a>!=B.<a class="code" href="classdcovector.html#a08e64e5ecb5121100787d0f5d48f0af4" title="vector size">l</a>){
<a name="l00266"></a>00266     ERROR_REPORT;
<a name="l00267"></a>00267     std::cerr &lt;&lt; <span class="stringliteral">&quot;These matrix and vector cannot be solved.&quot;</span> &lt;&lt; std::endl
<a name="l00268"></a>00268               &lt;&lt; <span class="stringliteral">&quot;Your input was (&quot;</span> &lt;&lt; <a class="code" href="classdgematrix.html#a3852aa1d0b876a74c1070e982cf388a2" title="matrix row size">m</a> &lt;&lt; <span class="stringliteral">&quot;x&quot;</span> &lt;&lt; <a class="code" href="classdgematrix.html#a60d26e3a499a235310e59a0aebbe0d0a" title="matrix column size">n</a> &lt;&lt; <span class="stringliteral">&quot;) and (&quot;</span> &lt;&lt; B.<a class="code" href="classdcovector.html#a08e64e5ecb5121100787d0f5d48f0af4" title="vector size">l</a> &lt;&lt; <span class="stringliteral">&quot;).&quot;</span> &lt;&lt; std::endl;
<a name="l00269"></a>00269     exit(1);
<a name="l00270"></a>00270   }
<a name="l00271"></a>00271 <span class="preprocessor">#endif//CPPL_DEBUG    </span>
<a name="l00272"></a>00272 <span class="preprocessor"></span>  
<a name="l00273"></a>00273   <span class="keywordflow">if</span>(<a class="code" href="classdgematrix.html#a3852aa1d0b876a74c1070e982cf388a2" title="matrix row size">m</a>&lt;<a class="code" href="classdgematrix.html#a60d26e3a499a235310e59a0aebbe0d0a" title="matrix column size">n</a>){ <span class="comment">//underdetermined</span>
<a name="l00274"></a>00274     <a class="code" href="classdcovector.html" title="Real Double-precision Column Vector Class.">dcovector</a> tmp(<a class="code" href="classdgematrix.html#a60d26e3a499a235310e59a0aebbe0d0a" title="matrix column size">n</a>);
<a name="l00275"></a>00275     <span class="keywordflow">for</span>(<span class="keywordtype">long</span> <a class="code" href="classdgematrix.html#a67f90160391aa6dbecc09b8eb6e23246">i</a>=0; <a class="code" href="classdgematrix.html#a67f90160391aa6dbecc09b8eb6e23246">i</a>&lt;B.<a class="code" href="classdcovector.html#a08e64e5ecb5121100787d0f5d48f0af4" title="vector size">l</a>; <a class="code" href="classdgematrix.html#a67f90160391aa6dbecc09b8eb6e23246">i</a>++){ tmp(<a class="code" href="classdgematrix.html#a67f90160391aa6dbecc09b8eb6e23246">i</a>)=B(<a class="code" href="classdgematrix.html#a67f90160391aa6dbecc09b8eb6e23246">i</a>); }
<a name="l00276"></a>00276     B.<a class="code" href="classdcovector.html#aa0a1e4ad0a058308b5c839609f05e973">clear</a>();
<a name="l00277"></a>00277     <a class="code" href="classdgematrix.html#a363aa32a0be4faacc1427b16aa59408c">swap</a>(B,tmp);
<a name="l00278"></a>00278   }
<a name="l00279"></a>00279   
<a name="l00280"></a>00280   S.<a class="code" href="classdcovector.html#a969d5aaa75cde3f7b1ddcfcdefddd7af">resize</a>(std::min(<a class="code" href="classdgematrix.html#a3852aa1d0b876a74c1070e982cf388a2" title="matrix row size">m</a>,<a class="code" href="classdgematrix.html#a60d26e3a499a235310e59a0aebbe0d0a" title="matrix column size">n</a>));
<a name="l00281"></a>00281   
<a name="l00282"></a>00282   <span class="keywordtype">long</span> NRHS(1), LDA(<a class="code" href="classdgematrix.html#a3852aa1d0b876a74c1070e982cf388a2" title="matrix row size">m</a>), LDB(B.<a class="code" href="classdcovector.html#a08e64e5ecb5121100787d0f5d48f0af4" title="vector size">l</a>),
<a name="l00283"></a>00283     LWORK(3*std::min(<a class="code" href="classdgematrix.html#a3852aa1d0b876a74c1070e982cf388a2" title="matrix row size">m</a>,<a class="code" href="classdgematrix.html#a60d26e3a499a235310e59a0aebbe0d0a" title="matrix column size">n</a>)+std::max(std::max(2*std::min(<a class="code" href="classdgematrix.html#a3852aa1d0b876a74c1070e982cf388a2" title="matrix row size">m</a>,<a class="code" href="classdgematrix.html#a60d26e3a499a235310e59a0aebbe0d0a" title="matrix column size">n</a>),std::max(<a class="code" href="classdgematrix.html#a3852aa1d0b876a74c1070e982cf388a2" title="matrix row size">m</a>,<a class="code" href="classdgematrix.html#a60d26e3a499a235310e59a0aebbe0d0a" title="matrix column size">n</a>)), NRHS)), INFO(1);
<a name="l00284"></a>00284   <span class="keywordtype">double</span> *WORK(<span class="keyword">new</span> <span class="keywordtype">double</span>[LWORK]);
<a name="l00285"></a>00285   dgelss_(<a class="code" href="classdgematrix.html#a3852aa1d0b876a74c1070e982cf388a2" title="matrix row size">m</a>, <a class="code" href="classdgematrix.html#a60d26e3a499a235310e59a0aebbe0d0a" title="matrix column size">n</a>, NRHS, <a class="code" href="classdgematrix.html#a006ce8c2f9add380bcb2816fdb13aa01" title="1D array to store matrix data">array</a>, LDA, B.<a class="code" href="classdcovector.html#a9b188bf160d8f7e24ea1d41deaceaf6d" title="1D array to store vector data">array</a>, LDB, S.<a class="code" href="classdcovector.html#a9b188bf160d8f7e24ea1d41deaceaf6d" title="1D array to store vector data">array</a>,
<a name="l00286"></a>00286           RCOND, RANK, WORK, LWORK, INFO);
<a name="l00287"></a>00287   <span class="keyword">delete</span> [] WORK;
<a name="l00288"></a>00288 
<a name="l00289"></a>00289   <span class="keywordflow">if</span>(INFO!=0){
<a name="l00290"></a>00290     WARNING_REPORT;
<a name="l00291"></a>00291     std::cerr &lt;&lt; <span class="stringliteral">&quot;Serious trouble happend. INFO = &quot;</span> &lt;&lt; INFO &lt;&lt; <span class="stringliteral">&quot;.&quot;</span> &lt;&lt; std::endl;
<a name="l00292"></a>00292   }
<a name="l00293"></a>00293   <span class="keywordflow">return</span> INFO;
<a name="l00294"></a>00294 }
<a name="l00295"></a>00295 
<a name="l00296"></a>00296 <span class="comment">//=============================================================================</span><span class="comment"></span>
<a name="l00297"></a>00297 <span class="comment">/*! calculate the least-squares-least-norm solution for overdetermined or </span>
<a name="l00298"></a>00298 <span class="comment">  underdetermined A*x=y using dgelss\n</span>
<a name="l00299"></a>00299 <span class="comment">*/</span>
<a name="l00300"></a><a class="code" href="classdgematrix.html#a99cdffa03cc68b2c1fd8a19490c2a443">00300</a> <span class="keyword">inline</span> <span class="keywordtype">long</span> <a class="code" href="classdgematrix.html#ab3abb914fecbe4b373b8bab0b4395230">dgematrix::dgelss</a>(<a class="code" href="classdgematrix.html" title="Real Double-precision General Dence Matrix Class.">dgematrix</a>&amp; B, <a class="code" href="classdcovector.html" title="Real Double-precision Column Vector Class.">dcovector</a>&amp; S, <span class="keywordtype">long</span>&amp; RANK,
<a name="l00301"></a>00301                               <span class="keyword">const</span> <span class="keywordtype">double</span> RCOND =-1. )
<a name="l00302"></a>00302 {VERBOSE_REPORT;
<a name="l00303"></a>00303 <span class="preprocessor">#ifdef  CPPL_DEBUG</span>
<a name="l00304"></a>00304 <span class="preprocessor"></span>  <span class="keywordflow">if</span>(<a class="code" href="classdgematrix.html#a3852aa1d0b876a74c1070e982cf388a2" title="matrix row size">m</a>!=B.<a class="code" href="classdgematrix.html#a3852aa1d0b876a74c1070e982cf388a2" title="matrix row size">m</a>){
<a name="l00305"></a>00305     ERROR_REPORT;
<a name="l00306"></a>00306     std::cerr &lt;&lt; <span class="stringliteral">&quot;These matrix and vector cannot be solved.&quot;</span> &lt;&lt; std::endl
<a name="l00307"></a>00307               &lt;&lt; <span class="stringliteral">&quot;Your input was (&quot;</span> &lt;&lt; <a class="code" href="classdgematrix.html#a3852aa1d0b876a74c1070e982cf388a2" title="matrix row size">m</a> &lt;&lt; <span class="stringliteral">&quot;x&quot;</span> &lt;&lt; <a class="code" href="classdgematrix.html#a60d26e3a499a235310e59a0aebbe0d0a" title="matrix column size">n</a> &lt;&lt; <span class="stringliteral">&quot;) and (&quot;</span> &lt;&lt; B.<a class="code" href="classdgematrix.html#a3852aa1d0b876a74c1070e982cf388a2" title="matrix row size">m</a> &lt;&lt; <span class="stringliteral">&quot;x&quot;</span> &lt;&lt; B.<a class="code" href="classdgematrix.html#a60d26e3a499a235310e59a0aebbe0d0a" title="matrix column size">n</a>  &lt;&lt; <span class="stringliteral">&quot;).&quot;</span> &lt;&lt; std::endl;
<a name="l00308"></a>00308     exit(1);
<a name="l00309"></a>00309   }
<a name="l00310"></a>00310 <span class="preprocessor">#endif//CPPL_DEBUG    </span>
<a name="l00311"></a>00311 <span class="preprocessor"></span>  
<a name="l00312"></a>00312   <span class="keywordflow">if</span>(<a class="code" href="classdgematrix.html#a3852aa1d0b876a74c1070e982cf388a2" title="matrix row size">m</a>&lt;<a class="code" href="classdgematrix.html#a60d26e3a499a235310e59a0aebbe0d0a" title="matrix column size">n</a>){ <span class="comment">//underdetermined</span>
<a name="l00313"></a>00313     <a class="code" href="classdgematrix.html" title="Real Double-precision General Dence Matrix Class.">dgematrix</a> tmp(<a class="code" href="classdgematrix.html#a60d26e3a499a235310e59a0aebbe0d0a" title="matrix column size">n</a>,B.<a class="code" href="classdgematrix.html#a60d26e3a499a235310e59a0aebbe0d0a" title="matrix column size">n</a>);
<a name="l00314"></a>00314     <span class="keywordflow">for</span>(<span class="keywordtype">long</span> <a class="code" href="classdgematrix.html#a67f90160391aa6dbecc09b8eb6e23246">i</a>=0; <a class="code" href="classdgematrix.html#a67f90160391aa6dbecc09b8eb6e23246">i</a>&lt;B.<a class="code" href="classdgematrix.html#a3852aa1d0b876a74c1070e982cf388a2" title="matrix row size">m</a>; <a class="code" href="classdgematrix.html#a67f90160391aa6dbecc09b8eb6e23246">i</a>++){
<a name="l00315"></a>00315       <span class="keywordflow">for</span>(<span class="keywordtype">long</span> j=0; j&lt;B.<a class="code" href="classdgematrix.html#a60d26e3a499a235310e59a0aebbe0d0a" title="matrix column size">n</a>; j++){
<a name="l00316"></a>00316         tmp(<a class="code" href="classdgematrix.html#a67f90160391aa6dbecc09b8eb6e23246">i</a>,j)=B(<a class="code" href="classdgematrix.html#a67f90160391aa6dbecc09b8eb6e23246">i</a>,j);
<a name="l00317"></a>00317       }
<a name="l00318"></a>00318     }
<a name="l00319"></a>00319     B.<a class="code" href="classdgematrix.html#aff82fa7ef089e6c565bfd4aaea4b9a27">clear</a>();
<a name="l00320"></a>00320     <a class="code" href="classdgematrix.html#a363aa32a0be4faacc1427b16aa59408c">swap</a>(B,tmp);
<a name="l00321"></a>00321   }
<a name="l00322"></a>00322   
<a name="l00323"></a>00323   S.<a class="code" href="classdcovector.html#a969d5aaa75cde3f7b1ddcfcdefddd7af">resize</a>(std::min(<a class="code" href="classdgematrix.html#a3852aa1d0b876a74c1070e982cf388a2" title="matrix row size">m</a>,<a class="code" href="classdgematrix.html#a60d26e3a499a235310e59a0aebbe0d0a" title="matrix column size">n</a>));
<a name="l00324"></a>00324   
<a name="l00325"></a>00325   <span class="keywordtype">long</span> NRHS(B.<a class="code" href="classdgematrix.html#a60d26e3a499a235310e59a0aebbe0d0a" title="matrix column size">n</a>), LDA(<a class="code" href="classdgematrix.html#a3852aa1d0b876a74c1070e982cf388a2" title="matrix row size">m</a>), LDB(B.<a class="code" href="classdgematrix.html#a3852aa1d0b876a74c1070e982cf388a2" title="matrix row size">m</a>),
<a name="l00326"></a>00326     LWORK(3*std::min(<a class="code" href="classdgematrix.html#a3852aa1d0b876a74c1070e982cf388a2" title="matrix row size">m</a>,<a class="code" href="classdgematrix.html#a60d26e3a499a235310e59a0aebbe0d0a" title="matrix column size">n</a>)+std::max(std::max(2*std::min(<a class="code" href="classdgematrix.html#a3852aa1d0b876a74c1070e982cf388a2" title="matrix row size">m</a>,<a class="code" href="classdgematrix.html#a60d26e3a499a235310e59a0aebbe0d0a" title="matrix column size">n</a>),std::max(<a class="code" href="classdgematrix.html#a3852aa1d0b876a74c1070e982cf388a2" title="matrix row size">m</a>,<a class="code" href="classdgematrix.html#a60d26e3a499a235310e59a0aebbe0d0a" title="matrix column size">n</a>)), NRHS)), INFO(1);
<a name="l00327"></a>00327   <span class="keywordtype">double</span> *WORK(<span class="keyword">new</span> <span class="keywordtype">double</span>[LWORK]);
<a name="l00328"></a>00328   dgelss_(<a class="code" href="classdgematrix.html#a3852aa1d0b876a74c1070e982cf388a2" title="matrix row size">m</a>, <a class="code" href="classdgematrix.html#a60d26e3a499a235310e59a0aebbe0d0a" title="matrix column size">n</a>, NRHS, <a class="code" href="classdgematrix.html#a006ce8c2f9add380bcb2816fdb13aa01" title="1D array to store matrix data">array</a>, LDA, B.<a class="code" href="classdgematrix.html#a006ce8c2f9add380bcb2816fdb13aa01" title="1D array to store matrix data">array</a>, LDB, S.<a class="code" href="classdcovector.html#a9b188bf160d8f7e24ea1d41deaceaf6d" title="1D array to store vector data">array</a>,
<a name="l00329"></a>00329           RCOND, RANK, WORK, LWORK, INFO);
<a name="l00330"></a>00330   <span class="keyword">delete</span> [] WORK;
<a name="l00331"></a>00331 
<a name="l00332"></a>00332   <span class="keywordflow">if</span>(INFO!=0){
<a name="l00333"></a>00333     WARNING_REPORT;
<a name="l00334"></a>00334     std::cerr &lt;&lt; <span class="stringliteral">&quot;Serious trouble happend. INFO = &quot;</span> &lt;&lt; INFO &lt;&lt; <span class="stringliteral">&quot;.&quot;</span> &lt;&lt; std::endl;
<a name="l00335"></a>00335   }
<a name="l00336"></a>00336   <span class="keywordflow">return</span> INFO;
<a name="l00337"></a>00337 }
<a name="l00338"></a>00338 <span class="comment"></span>
<a name="l00339"></a>00339 <span class="comment">///////////////////////////////////////////////////////////////////////////////</span>
<a name="l00340"></a>00340 <span class="comment">///////////////////////////////////////////////////////////////////////////////</span>
<a name="l00341"></a>00341 <span class="comment">///////////////////////////////////////////////////////////////////////////////</span>
<a name="l00342"></a>00342 <span class="comment"></span>
<a name="l00343"></a>00343 <span class="comment">//=============================================================================</span><span class="comment"></span>
<a name="l00344"></a>00344 <span class="comment">/*! calculate eigenvalues\n</span>
<a name="l00345"></a>00345 <span class="comment">  All of the arguments need not to be initialized. </span>
<a name="l00346"></a>00346 <span class="comment">  wr and wi are overwitten and become </span>
<a name="l00347"></a>00347 <span class="comment">  real and imaginary part of eigenvalues, respectively. </span>
<a name="l00348"></a>00348 <span class="comment">  This matrix is also overwritten. </span>
<a name="l00349"></a>00349 <span class="comment">*/</span>
<a name="l00350"></a><a class="code" href="classdgematrix.html#ab4c6d8212ee192595eca196890defb81">00350</a> <span class="keyword">inline</span> <span class="keywordtype">long</span> <a class="code" href="classdgematrix.html#ab4c6d8212ee192595eca196890defb81">dgematrix::dgeev</a>(std::vector&lt;double&gt;&amp; wr, std::vector&lt;double&gt;&amp; wi)
<a name="l00351"></a>00351 {VERBOSE_REPORT;
<a name="l00352"></a>00352 <span class="preprocessor">#ifdef  CPPL_DEBUG</span>
<a name="l00353"></a>00353 <span class="preprocessor"></span>  <span class="keywordflow">if</span>(<a class="code" href="classdgematrix.html#a3852aa1d0b876a74c1070e982cf388a2" title="matrix row size">m</a>!=<a class="code" href="classdgematrix.html#a60d26e3a499a235310e59a0aebbe0d0a" title="matrix column size">n</a>){
<a name="l00354"></a>00354     ERROR_REPORT;
<a name="l00355"></a>00355     std::cerr &lt;&lt; <span class="stringliteral">&quot;This matrix is not a square matrix.&quot;</span> &lt;&lt; std::endl
<a name="l00356"></a>00356               &lt;&lt; <span class="stringliteral">&quot;This matrix is (&quot;</span> &lt;&lt; <a class="code" href="classdgematrix.html#a3852aa1d0b876a74c1070e982cf388a2" title="matrix row size">m</a> &lt;&lt; <span class="stringliteral">&quot;x&quot;</span> &lt;&lt; <a class="code" href="classdgematrix.html#a60d26e3a499a235310e59a0aebbe0d0a" title="matrix column size">n</a> &lt;&lt; <span class="stringliteral">&quot;).&quot;</span> &lt;&lt; std::endl;
<a name="l00357"></a>00357     exit(1);
<a name="l00358"></a>00358   }
<a name="l00359"></a>00359 <span class="preprocessor">#endif//CPPL_DEBUG</span>
<a name="l00360"></a>00360 <span class="preprocessor"></span>  
<a name="l00361"></a>00361   wr.resize(<a class="code" href="classdgematrix.html#a60d26e3a499a235310e59a0aebbe0d0a" title="matrix column size">n</a>); wi.resize(<a class="code" href="classdgematrix.html#a60d26e3a499a235310e59a0aebbe0d0a" title="matrix column size">n</a>);
<a name="l00362"></a>00362   <span class="keywordtype">char</span> JOBVL(<span class="charliteral">&#39;n&#39;</span>), JOBVR(<span class="charliteral">&#39;n&#39;</span>);
<a name="l00363"></a>00363   <span class="keywordtype">long</span> LDA(<a class="code" href="classdgematrix.html#a60d26e3a499a235310e59a0aebbe0d0a" title="matrix column size">n</a>), LDVL(1), LDVR(1), LWORK(3*<a class="code" href="classdgematrix.html#a60d26e3a499a235310e59a0aebbe0d0a" title="matrix column size">n</a>), INFO(1);
<a name="l00364"></a>00364   <span class="keywordtype">double</span> *VL(NULL), *VR(NULL), *WORK(<span class="keyword">new</span> <span class="keywordtype">double</span>[LWORK]);
<a name="l00365"></a>00365   dgeev_(JOBVL, JOBVR, <a class="code" href="classdgematrix.html#a60d26e3a499a235310e59a0aebbe0d0a" title="matrix column size">n</a>, <a class="code" href="classdgematrix.html#a006ce8c2f9add380bcb2816fdb13aa01" title="1D array to store matrix data">array</a>, LDA, &amp;wr[0], &amp;wi[0], 
<a name="l00366"></a>00366          VL, LDVL, VR, LDVR, WORK, LWORK, INFO);
<a name="l00367"></a>00367   <span class="keyword">delete</span> [] WORK; <span class="keyword">delete</span> [] VL; <span class="keyword">delete</span> [] VL;
<a name="l00368"></a>00368   
<a name="l00369"></a>00369   <span class="keywordflow">if</span>(INFO!=0){
<a name="l00370"></a>00370     WARNING_REPORT;
<a name="l00371"></a>00371     std::cerr &lt;&lt; <span class="stringliteral">&quot;Serious trouble happend. INFO = &quot;</span> &lt;&lt; INFO &lt;&lt; <span class="stringliteral">&quot;.&quot;</span> &lt;&lt; std::endl;
<a name="l00372"></a>00372   }
<a name="l00373"></a>00373   <span class="keywordflow">return</span> INFO;
<a name="l00374"></a>00374 }
<a name="l00375"></a>00375 
<a name="l00376"></a>00376 <span class="comment">//=============================================================================</span><span class="comment"></span>
<a name="l00377"></a>00377 <span class="comment">/*! calculate eigenvalues\n</span>
<a name="l00378"></a>00378 <span class="comment">  All of the arguments need not to be initialized. </span>
<a name="l00379"></a>00379 <span class="comment">  w are overwitten and become eigenvalues, respectively.</span>
<a name="l00380"></a>00380 <span class="comment">  This matrix is also overwritten. </span>
<a name="l00381"></a>00381 <span class="comment">*/</span>
<a name="l00382"></a><a class="code" href="classdgematrix.html#ac00b17d31d60895c2786d030b816a930">00382</a> <span class="keyword">inline</span> <span class="keywordtype">long</span> <a class="code" href="classdgematrix.html#ab4c6d8212ee192595eca196890defb81">dgematrix::dgeev</a>(<a class="code" href="classzcovector.html" title="Complex Double-precision Column Vector Class.">zcovector</a>&amp; w)
<a name="l00383"></a>00383 {VERBOSE_REPORT;<span class="comment"></span>
<a name="l00384"></a>00384 <span class="comment">  //// call dgeev ////</span>
<a name="l00385"></a>00385 <span class="comment"></span>  std::vector&lt;double&gt; wr, wi;
<a name="l00386"></a>00386   <span class="keywordtype">long</span> INFO =<a class="code" href="classdgematrix.html#ab4c6d8212ee192595eca196890defb81">dgeev</a>(wr,wi);
<a name="l00387"></a>00387   <span class="comment"></span>
<a name="l00388"></a>00388 <span class="comment">  //// assign ////</span>
<a name="l00389"></a>00389 <span class="comment"></span>  w.<a class="code" href="classzcovector.html#ae1cbcca07337389700b0837842094792">resize</a>(<a class="code" href="classdgematrix.html#a60d26e3a499a235310e59a0aebbe0d0a" title="matrix column size">n</a>);
<a name="l00390"></a>00390   <span class="keywordflow">for</span>(<span class="keywordtype">long</span> <a class="code" href="classdgematrix.html#a67f90160391aa6dbecc09b8eb6e23246">i</a>=0; <a class="code" href="classdgematrix.html#a67f90160391aa6dbecc09b8eb6e23246">i</a>&lt;<a class="code" href="classdgematrix.html#a60d26e3a499a235310e59a0aebbe0d0a" title="matrix column size">n</a>; <a class="code" href="classdgematrix.html#a67f90160391aa6dbecc09b8eb6e23246">i</a>++){
<a name="l00391"></a>00391     w(<a class="code" href="classdgematrix.html#a67f90160391aa6dbecc09b8eb6e23246">i</a>) =comple(wr[<a class="code" href="classdgematrix.html#a67f90160391aa6dbecc09b8eb6e23246">i</a>],wi[i]);
<a name="l00392"></a>00392   }
<a name="l00393"></a>00393   
<a name="l00394"></a>00394   <span class="keywordflow">return</span> INFO;
<a name="l00395"></a>00395 }
<a name="l00396"></a>00396 
<a name="l00397"></a>00397 <span class="comment">//=============================================================================</span><span class="comment"></span>
<a name="l00398"></a>00398 <span class="comment">/*! calculate right eigenvalues and right eigenvectors\n</span>
<a name="l00399"></a>00399 <span class="comment">  All of the arguments need not to be initialized. </span>
<a name="l00400"></a>00400 <span class="comment">  wr, wi, vrr, vri are overwitten and become </span>
<a name="l00401"></a>00401 <span class="comment">  real and imaginary part of right eigenvalues and right eigenvectors, </span>
<a name="l00402"></a>00402 <span class="comment">  respectively. </span>
<a name="l00403"></a>00403 <span class="comment">  This matrix is also overwritten. </span>
<a name="l00404"></a>00404 <span class="comment">*/</span>
<a name="l00405"></a>00405 <span class="keyword">inline</span> <span class="keywordtype">long</span> <a class="code" href="classdgematrix.html#ab4c6d8212ee192595eca196890defb81">dgematrix::dgeev</a>
<a name="l00406"></a><a class="code" href="classdgematrix.html#a82152e07a32a87a430dccdb9f93f3ac0">00406</a> (
<a name="l00407"></a>00407  std::vector&lt;double&gt;&amp; wr,
<a name="l00408"></a>00408  std::vector&lt;double&gt;&amp; wi, 
<a name="l00409"></a>00409  std::vector&lt;dcovector&gt;&amp; vrr,
<a name="l00410"></a>00410  std::vector&lt;dcovector&gt;&amp; vri
<a name="l00411"></a>00411  )
<a name="l00412"></a>00412 {VERBOSE_REPORT;
<a name="l00413"></a>00413 <span class="preprocessor">#ifdef  CPPL_DEBUG</span>
<a name="l00414"></a>00414 <span class="preprocessor"></span>  <span class="keywordflow">if</span>(m!=n){
<a name="l00415"></a>00415     ERROR_REPORT;
<a name="l00416"></a>00416     std::cerr &lt;&lt; <span class="stringliteral">&quot;This matrix is not a square matrix.&quot;</span> &lt;&lt; std::endl
<a name="l00417"></a>00417               &lt;&lt; <span class="stringliteral">&quot;This matrix is (&quot;</span> &lt;&lt; m &lt;&lt; <span class="stringliteral">&quot;x&quot;</span> &lt;&lt; n &lt;&lt; <span class="stringliteral">&quot;).&quot;</span> &lt;&lt; std::endl;
<a name="l00418"></a>00418     exit(1);
<a name="l00419"></a>00419   }
<a name="l00420"></a>00420 <span class="preprocessor">#endif//CPPL_DEBUG</span>
<a name="l00421"></a>00421 <span class="preprocessor"></span>  
<a name="l00422"></a>00422   wr.resize(n); wi.resize(n); vrr.resize(n); vri.resize(n);
<a name="l00423"></a>00423   <span class="keywordflow">for</span>(<span class="keywordtype">long</span> <a class="code" href="__dgbmatrix-calc_8hpp.html#af7c79b3740ba9c1d547f55ac94bc1258">i</a>=0; <a class="code" href="__dgbmatrix-calc_8hpp.html#af7c79b3740ba9c1d547f55ac94bc1258">i</a>&lt;n; <a class="code" href="__dgbmatrix-calc_8hpp.html#af7c79b3740ba9c1d547f55ac94bc1258">i</a>++){ vrr[<a class="code" href="__dgbmatrix-calc_8hpp.html#af7c79b3740ba9c1d547f55ac94bc1258">i</a>].resize(n); vri[<a class="code" href="__dgbmatrix-calc_8hpp.html#af7c79b3740ba9c1d547f55ac94bc1258">i</a>].resize(n); }
<a name="l00424"></a>00424   <a class="code" href="classdgematrix.html" title="Real Double-precision General Dence Matrix Class.">dgematrix</a> VR(n,n);
<a name="l00425"></a>00425   <span class="keywordtype">char</span> JOBVL(<span class="charliteral">&#39;n&#39;</span>), JOBVR(<span class="charliteral">&#39;V&#39;</span>);
<a name="l00426"></a>00426   <span class="keywordtype">long</span> LDA(n), LDVL(1), LDVR(n), LWORK(4*n), INFO(1);
<a name="l00427"></a>00427   <span class="keywordtype">double</span> *VL(NULL), *WORK(<span class="keyword">new</span> <span class="keywordtype">double</span>[LWORK]);
<a name="l00428"></a>00428   dgeev_(JOBVL, JOBVR, n, array, LDA, &amp;wr[0], &amp;wi[0], 
<a name="l00429"></a>00429          VL, LDVL, VR.<a class="code" href="classdgematrix.html#a006ce8c2f9add380bcb2816fdb13aa01" title="1D array to store matrix data">array</a>, LDVR, WORK, LWORK, INFO);
<a name="l00430"></a>00430   <span class="keyword">delete</span> [] WORK; <span class="keyword">delete</span> [] VL;
<a name="l00431"></a>00431   <span class="comment"></span>
<a name="l00432"></a>00432 <span class="comment">  //// forming ////</span>
<a name="l00433"></a>00433 <span class="comment"></span>  <span class="keywordflow">for</span>(<span class="keywordtype">long</span> j=0; j&lt;n; j++){
<a name="l00434"></a>00434     <span class="keywordflow">if</span>(fabs(wi[j])&lt;DBL_MIN){
<a name="l00435"></a>00435       <span class="keywordflow">for</span>(<span class="keywordtype">long</span> <a class="code" href="__dgbmatrix-calc_8hpp.html#af7c79b3740ba9c1d547f55ac94bc1258">i</a>=0; <a class="code" href="__dgbmatrix-calc_8hpp.html#af7c79b3740ba9c1d547f55ac94bc1258">i</a>&lt;n; <a class="code" href="__dgbmatrix-calc_8hpp.html#af7c79b3740ba9c1d547f55ac94bc1258">i</a>++){
<a name="l00436"></a>00436         vrr[j](<a class="code" href="__dgbmatrix-calc_8hpp.html#af7c79b3740ba9c1d547f55ac94bc1258">i</a>) = VR(<a class="code" href="__dgbmatrix-calc_8hpp.html#af7c79b3740ba9c1d547f55ac94bc1258">i</a>,j);
<a name="l00437"></a>00437         vri[j](<a class="code" href="__dgbmatrix-calc_8hpp.html#af7c79b3740ba9c1d547f55ac94bc1258">i</a>) = 0.0;
<a name="l00438"></a>00438       }
<a name="l00439"></a>00439     }
<a name="l00440"></a>00440     <span class="keywordflow">else</span>{
<a name="l00441"></a>00441       <span class="keywordflow">for</span>(<span class="keywordtype">long</span> <a class="code" href="__dgbmatrix-calc_8hpp.html#af7c79b3740ba9c1d547f55ac94bc1258">i</a>=0; <a class="code" href="__dgbmatrix-calc_8hpp.html#af7c79b3740ba9c1d547f55ac94bc1258">i</a>&lt;n; <a class="code" href="__dgbmatrix-calc_8hpp.html#af7c79b3740ba9c1d547f55ac94bc1258">i</a>++){
<a name="l00442"></a>00442         vrr[j](<a class="code" href="__dgbmatrix-calc_8hpp.html#af7c79b3740ba9c1d547f55ac94bc1258">i</a>)   = VR(<a class="code" href="__dgbmatrix-calc_8hpp.html#af7c79b3740ba9c1d547f55ac94bc1258">i</a>,j);
<a name="l00443"></a>00443         vri[j](<a class="code" href="__dgbmatrix-calc_8hpp.html#af7c79b3740ba9c1d547f55ac94bc1258">i</a>)   = VR(<a class="code" href="__dgbmatrix-calc_8hpp.html#af7c79b3740ba9c1d547f55ac94bc1258">i</a>,j+1);
<a name="l00444"></a>00444         vrr[j+1](<a class="code" href="__dgbmatrix-calc_8hpp.html#af7c79b3740ba9c1d547f55ac94bc1258">i</a>) = VR(<a class="code" href="__dgbmatrix-calc_8hpp.html#af7c79b3740ba9c1d547f55ac94bc1258">i</a>,j);
<a name="l00445"></a>00445         vri[j+1](<a class="code" href="__dgbmatrix-calc_8hpp.html#af7c79b3740ba9c1d547f55ac94bc1258">i</a>) =-VR(<a class="code" href="__dgbmatrix-calc_8hpp.html#af7c79b3740ba9c1d547f55ac94bc1258">i</a>,j+1);
<a name="l00446"></a>00446       }
<a name="l00447"></a>00447       j++;
<a name="l00448"></a>00448     }
<a name="l00449"></a>00449   }
<a name="l00450"></a>00450   
<a name="l00451"></a>00451   <span class="keywordflow">if</span>(INFO!=0){
<a name="l00452"></a>00452     WARNING_REPORT;
<a name="l00453"></a>00453     std::cerr &lt;&lt; <span class="stringliteral">&quot;Serious trouble happend. INFO = &quot;</span> &lt;&lt; INFO &lt;&lt; <span class="stringliteral">&quot;.&quot;</span> &lt;&lt; std::endl;
<a name="l00454"></a>00454   }
<a name="l00455"></a>00455   <span class="keywordflow">return</span> INFO;
<a name="l00456"></a>00456 }
<a name="l00457"></a>00457 
<a name="l00458"></a>00458 <span class="comment">//=============================================================================</span><span class="comment"></span>
<a name="l00459"></a>00459 <span class="comment">/*! calculate left eigenvalues and left eigenvectors\n</span>
<a name="l00460"></a>00460 <span class="comment">  All of the arguments need not to be initialized. </span>
<a name="l00461"></a>00461 <span class="comment">  wr, wi, vrr, vri are overwitten and become </span>
<a name="l00462"></a>00462 <span class="comment">  real and imaginary part of left eigenvalues and left eigenvectors, </span>
<a name="l00463"></a>00463 <span class="comment">  respectively. </span>
<a name="l00464"></a>00464 <span class="comment">  This matrix is also overwritten. </span>
<a name="l00465"></a>00465 <span class="comment">*/</span>
<a name="l00466"></a><a class="code" href="classdgematrix.html#adbabdbe158f8c57c2ba95515f0590105">00466</a> <span class="keyword">inline</span> <span class="keywordtype">long</span> <a class="code" href="classdgematrix.html#ab4c6d8212ee192595eca196890defb81">dgematrix::dgeev</a>(std::vector&lt;double&gt;&amp; wr, std::vector&lt;double&gt;&amp; wi, 
<a name="l00467"></a>00467                              std::vector&lt;drovector&gt;&amp; vlr, 
<a name="l00468"></a>00468                              std::vector&lt;drovector&gt;&amp; vli)
<a name="l00469"></a>00469 {VERBOSE_REPORT;
<a name="l00470"></a>00470 <span class="preprocessor">#ifdef  CPPL_DEBUG</span>
<a name="l00471"></a>00471 <span class="preprocessor"></span>  <span class="keywordflow">if</span>(<a class="code" href="classdgematrix.html#a3852aa1d0b876a74c1070e982cf388a2" title="matrix row size">m</a>!=<a class="code" href="classdgematrix.html#a60d26e3a499a235310e59a0aebbe0d0a" title="matrix column size">n</a>){
<a name="l00472"></a>00472     ERROR_REPORT;
<a name="l00473"></a>00473     std::cerr &lt;&lt; <span class="stringliteral">&quot;This matrix is not a square matrix.&quot;</span> &lt;&lt; std::endl
<a name="l00474"></a>00474               &lt;&lt; <span class="stringliteral">&quot;This matrix is (&quot;</span> &lt;&lt; <a class="code" href="classdgematrix.html#a3852aa1d0b876a74c1070e982cf388a2" title="matrix row size">m</a> &lt;&lt; <span class="stringliteral">&quot;x&quot;</span> &lt;&lt; <a class="code" href="classdgematrix.html#a60d26e3a499a235310e59a0aebbe0d0a" title="matrix column size">n</a> &lt;&lt; <span class="stringliteral">&quot;).&quot;</span> &lt;&lt; std::endl;
<a name="l00475"></a>00475     exit(1);
<a name="l00476"></a>00476   }
<a name="l00477"></a>00477 <span class="preprocessor">#endif//CPPL_DEBUG</span>
<a name="l00478"></a>00478 <span class="preprocessor"></span>  
<a name="l00479"></a>00479   wr.resize(<a class="code" href="classdgematrix.html#a60d26e3a499a235310e59a0aebbe0d0a" title="matrix column size">n</a>); wi.resize(<a class="code" href="classdgematrix.html#a60d26e3a499a235310e59a0aebbe0d0a" title="matrix column size">n</a>); vlr.resize(<a class="code" href="classdgematrix.html#a60d26e3a499a235310e59a0aebbe0d0a" title="matrix column size">n</a>); vli.resize(<a class="code" href="classdgematrix.html#a60d26e3a499a235310e59a0aebbe0d0a" title="matrix column size">n</a>);
<a name="l00480"></a>00480   <span class="keywordflow">for</span>(<span class="keywordtype">long</span> <a class="code" href="classdgematrix.html#a67f90160391aa6dbecc09b8eb6e23246">i</a>=0; <a class="code" href="classdgematrix.html#a67f90160391aa6dbecc09b8eb6e23246">i</a>&lt;<a class="code" href="classdgematrix.html#a60d26e3a499a235310e59a0aebbe0d0a" title="matrix column size">n</a>; <a class="code" href="classdgematrix.html#a67f90160391aa6dbecc09b8eb6e23246">i</a>++){ vlr[<a class="code" href="classdgematrix.html#a67f90160391aa6dbecc09b8eb6e23246">i</a>].resize(n); vli[<a class="code" href="classdgematrix.html#a67f90160391aa6dbecc09b8eb6e23246">i</a>].resize(n); }
<a name="l00481"></a>00481   <a class="code" href="classdgematrix.html" title="Real Double-precision General Dence Matrix Class.">dgematrix</a> VL(n,n);
<a name="l00482"></a>00482   <span class="keywordtype">char</span> JOBVL(<span class="charliteral">&#39;V&#39;</span>), JOBVR(<span class="charliteral">&#39;n&#39;</span>);
<a name="l00483"></a>00483   <span class="keywordtype">long</span> LDA(n), LDVL(n), LDVR(1), LWORK(4*n), INFO(1);
<a name="l00484"></a>00484   <span class="keywordtype">double</span> *VR(NULL), *WORK(<span class="keyword">new</span> <span class="keywordtype">double</span>[LWORK]);
<a name="l00485"></a>00485   dgeev_(JOBVL, JOBVR, n, <a class="code" href="classdgematrix.html#a006ce8c2f9add380bcb2816fdb13aa01" title="1D array to store matrix data">array</a>, LDA, &amp;wr[0], &amp;wi[0], 
<a name="l00486"></a>00486          VL.<a class="code" href="classdgematrix.html#a006ce8c2f9add380bcb2816fdb13aa01" title="1D array to store matrix data">array</a>, LDVL, VR, LDVR, WORK, LWORK, INFO);
<a name="l00487"></a>00487   <span class="keyword">delete</span> [] WORK; <span class="keyword">delete</span> [] VR;
<a name="l00488"></a>00488 <span class="comment"></span>
<a name="l00489"></a>00489 <span class="comment">  //// forming ////</span>
<a name="l00490"></a>00490 <span class="comment"></span>  <span class="keywordflow">for</span>(<span class="keywordtype">long</span> j=0; j&lt;<a class="code" href="classdgematrix.html#a60d26e3a499a235310e59a0aebbe0d0a" title="matrix column size">n</a>; j++){
<a name="l00491"></a>00491     <span class="keywordflow">if</span>(fabs(wi[j])&lt;DBL_MIN){
<a name="l00492"></a>00492       <span class="keywordflow">for</span>(<span class="keywordtype">long</span> <a class="code" href="classdgematrix.html#a67f90160391aa6dbecc09b8eb6e23246">i</a>=0; <a class="code" href="classdgematrix.html#a67f90160391aa6dbecc09b8eb6e23246">i</a>&lt;<a class="code" href="classdgematrix.html#a60d26e3a499a235310e59a0aebbe0d0a" title="matrix column size">n</a>; <a class="code" href="classdgematrix.html#a67f90160391aa6dbecc09b8eb6e23246">i</a>++){
<a name="l00493"></a>00493         vlr[j](<a class="code" href="classdgematrix.html#a67f90160391aa6dbecc09b8eb6e23246">i</a>) = VL(<a class="code" href="classdgematrix.html#a67f90160391aa6dbecc09b8eb6e23246">i</a>,j);
<a name="l00494"></a>00494         vli[j](<a class="code" href="classdgematrix.html#a67f90160391aa6dbecc09b8eb6e23246">i</a>) = 0.0;
<a name="l00495"></a>00495       }
<a name="l00496"></a>00496     }
<a name="l00497"></a>00497     <span class="keywordflow">else</span>{
<a name="l00498"></a>00498       <span class="keywordflow">for</span>(<span class="keywordtype">long</span> <a class="code" href="classdgematrix.html#a67f90160391aa6dbecc09b8eb6e23246">i</a>=0; <a class="code" href="classdgematrix.html#a67f90160391aa6dbecc09b8eb6e23246">i</a>&lt;<a class="code" href="classdgematrix.html#a60d26e3a499a235310e59a0aebbe0d0a" title="matrix column size">n</a>; <a class="code" href="classdgematrix.html#a67f90160391aa6dbecc09b8eb6e23246">i</a>++){
<a name="l00499"></a>00499         vlr[j](<a class="code" href="classdgematrix.html#a67f90160391aa6dbecc09b8eb6e23246">i</a>)   = VL(<a class="code" href="classdgematrix.html#a67f90160391aa6dbecc09b8eb6e23246">i</a>,j);
<a name="l00500"></a>00500         vli[j](<a class="code" href="classdgematrix.html#a67f90160391aa6dbecc09b8eb6e23246">i</a>)   =-VL(<a class="code" href="classdgematrix.html#a67f90160391aa6dbecc09b8eb6e23246">i</a>,j+1);
<a name="l00501"></a>00501         vlr[j+1](<a class="code" href="classdgematrix.html#a67f90160391aa6dbecc09b8eb6e23246">i</a>) = VL(<a class="code" href="classdgematrix.html#a67f90160391aa6dbecc09b8eb6e23246">i</a>,j);
<a name="l00502"></a>00502         vli[j+1](<a class="code" href="classdgematrix.html#a67f90160391aa6dbecc09b8eb6e23246">i</a>) = VL(<a class="code" href="classdgematrix.html#a67f90160391aa6dbecc09b8eb6e23246">i</a>,j+1);
<a name="l00503"></a>00503       }
<a name="l00504"></a>00504       j++;
<a name="l00505"></a>00505     }
<a name="l00506"></a>00506   }
<a name="l00507"></a>00507   
<a name="l00508"></a>00508 
<a name="l00509"></a>00509   <span class="keywordflow">if</span>(INFO!=0){
<a name="l00510"></a>00510     WARNING_REPORT;
<a name="l00511"></a>00511     std::cerr &lt;&lt; <span class="stringliteral">&quot;Serious trouble happend. INFO = &quot;</span> &lt;&lt; INFO &lt;&lt; <span class="stringliteral">&quot;.&quot;</span> &lt;&lt; std::endl;
<a name="l00512"></a>00512   }
<a name="l00513"></a>00513   <span class="keywordflow">return</span> INFO;
<a name="l00514"></a>00514 }
<a name="l00515"></a>00515 
<a name="l00516"></a>00516 
<a name="l00517"></a>00517 <span class="comment"></span>
<a name="l00518"></a>00518 <span class="comment">///////////////////////////////////////////////////////////////////////////////</span>
<a name="l00519"></a>00519 <span class="comment">///////////////////////////////////////////////////////////////////////////////</span>
<a name="l00520"></a>00520 <span class="comment">///////////////////////////////////////////////////////////////////////////////</span>
<a name="l00521"></a>00521 <span class="comment"></span>
<a name="l00522"></a>00522 <span class="comment">//=============================================================================</span><span class="comment"></span>
<a name="l00523"></a>00523 <span class="comment">/*! calculate generalized eigenvalues\n</span>
<a name="l00524"></a>00524 <span class="comment">  All of the arguments don&#39;t need to be initialized. </span>
<a name="l00525"></a>00525 <span class="comment">  wr and wi are overwitten and become </span>
<a name="l00526"></a>00526 <span class="comment">  real and imaginary part of generalized eigenvalues, respectively. </span>
<a name="l00527"></a>00527 <span class="comment">  This matrix and matB are also overwritten. </span>
<a name="l00528"></a>00528 <span class="comment">*/</span>
<a name="l00529"></a><a class="code" href="classdgematrix.html#a3374d63da801c823055f52dc7f2f55d9">00529</a> <span class="keyword">inline</span> <span class="keywordtype">long</span> <a class="code" href="classdgematrix.html#a3374d63da801c823055f52dc7f2f55d9">dgematrix::dggev</a>(<a class="code" href="classdgematrix.html" title="Real Double-precision General Dence Matrix Class.">dgematrix</a>&amp; matB,
<a name="l00530"></a>00530                              std::vector&lt;double&gt;&amp; wr, std::vector&lt;double&gt;&amp; wi)
<a name="l00531"></a>00531 {VERBOSE_REPORT;
<a name="l00532"></a>00532 <span class="preprocessor">#ifdef  CPPL_DEBUG</span>
<a name="l00533"></a>00533 <span class="preprocessor"></span>  <span class="keywordflow">if</span>(<a class="code" href="classdgematrix.html#a3852aa1d0b876a74c1070e982cf388a2" title="matrix row size">m</a>!=<a class="code" href="classdgematrix.html#a60d26e3a499a235310e59a0aebbe0d0a" title="matrix column size">n</a>){
<a name="l00534"></a>00534     ERROR_REPORT;
<a name="l00535"></a>00535     std::cerr &lt;&lt; <span class="stringliteral">&quot;This matrix is not a square matrix.&quot;</span> &lt;&lt; std::endl
<a name="l00536"></a>00536               &lt;&lt; <span class="stringliteral">&quot;This matrix is (&quot;</span> &lt;&lt; <a class="code" href="classdgematrix.html#a3852aa1d0b876a74c1070e982cf388a2" title="matrix row size">m</a> &lt;&lt; <span class="stringliteral">&quot;x&quot;</span> &lt;&lt; <a class="code" href="classdgematrix.html#a60d26e3a499a235310e59a0aebbe0d0a" title="matrix column size">n</a> &lt;&lt; <span class="stringliteral">&quot;).&quot;</span> &lt;&lt; std::endl;
<a name="l00537"></a>00537     exit(1);
<a name="l00538"></a>00538   }
<a name="l00539"></a>00539   <span class="keywordflow">if</span>(matB.<a class="code" href="classdgematrix.html#a3852aa1d0b876a74c1070e982cf388a2" title="matrix row size">m</a>!=<a class="code" href="classdgematrix.html#a60d26e3a499a235310e59a0aebbe0d0a" title="matrix column size">n</a> || matB.<a class="code" href="classdgematrix.html#a60d26e3a499a235310e59a0aebbe0d0a" title="matrix column size">n</a>!=<a class="code" href="classdgematrix.html#a60d26e3a499a235310e59a0aebbe0d0a" title="matrix column size">n</a>){
<a name="l00540"></a>00540     ERROR_REPORT;
<a name="l00541"></a>00541     std::cerr &lt;&lt; <span class="stringliteral">&quot;The matrix B is not a square matrix having the same size as \&quot;this\&quot; matrix.&quot;</span> &lt;&lt; std::endl
<a name="l00542"></a>00542               &lt;&lt; <span class="stringliteral">&quot;The B matrix is (&quot;</span> &lt;&lt; matB.<a class="code" href="classdgematrix.html#a3852aa1d0b876a74c1070e982cf388a2" title="matrix row size">m</a> &lt;&lt; <span class="stringliteral">&quot;x&quot;</span> &lt;&lt; matB.<a class="code" href="classdgematrix.html#a60d26e3a499a235310e59a0aebbe0d0a" title="matrix column size">n</a> &lt;&lt; <span class="stringliteral">&quot;).&quot;</span> &lt;&lt; std::endl;
<a name="l00543"></a>00543     exit(1);
<a name="l00544"></a>00544   }
<a name="l00545"></a>00545 <span class="preprocessor">#endif//CPPL_DEBUG</span>
<a name="l00546"></a>00546 <span class="preprocessor"></span>  
<a name="l00547"></a>00547   wr.resize(<a class="code" href="classdgematrix.html#a60d26e3a499a235310e59a0aebbe0d0a" title="matrix column size">n</a>); wi.resize(<a class="code" href="classdgematrix.html#a60d26e3a499a235310e59a0aebbe0d0a" title="matrix column size">n</a>);
<a name="l00548"></a>00548   <span class="keywordtype">char</span> JOBVL(<span class="charliteral">&#39;n&#39;</span>), JOBVR(<span class="charliteral">&#39;n&#39;</span>);
<a name="l00549"></a>00549   <span class="keywordtype">long</span> LDA(<a class="code" href="classdgematrix.html#a60d26e3a499a235310e59a0aebbe0d0a" title="matrix column size">n</a>), LDB(<a class="code" href="classdgematrix.html#a60d26e3a499a235310e59a0aebbe0d0a" title="matrix column size">n</a>), LDVL(1), LDVR(1), LWORK(8*<a class="code" href="classdgematrix.html#a60d26e3a499a235310e59a0aebbe0d0a" title="matrix column size">n</a>), INFO(1);
<a name="l00550"></a>00550   <span class="keywordtype">double</span> *BETA(<span class="keyword">new</span> <span class="keywordtype">double</span>[<a class="code" href="classdgematrix.html#a60d26e3a499a235310e59a0aebbe0d0a" title="matrix column size">n</a>]), *VL(NULL), *VR(NULL), 
<a name="l00551"></a>00551     *WORK(<span class="keyword">new</span> <span class="keywordtype">double</span>[LWORK]);
<a name="l00552"></a>00552   dggev_(JOBVL, JOBVR, n, <a class="code" href="classdgematrix.html#a006ce8c2f9add380bcb2816fdb13aa01" title="1D array to store matrix data">array</a>, LDA, matB.<a class="code" href="classdgematrix.html#a006ce8c2f9add380bcb2816fdb13aa01" title="1D array to store matrix data">array</a>, LDB, &amp;wr[0], &amp;wi[0], 
<a name="l00553"></a>00553          BETA, VL, LDVL, VR, LDVR, WORK, LWORK, INFO);
<a name="l00554"></a>00554   <span class="keyword">delete</span> [] WORK; <span class="keyword">delete</span> [] VL; <span class="keyword">delete</span> [] VR;
<a name="l00555"></a>00555   <span class="comment"></span>
<a name="l00556"></a>00556 <span class="comment">  //// reforming ////</span>
<a name="l00557"></a>00557 <span class="comment"></span>  <span class="keywordflow">for</span>(<span class="keywordtype">long</span> <a class="code" href="classdgematrix.html#a67f90160391aa6dbecc09b8eb6e23246">i</a>=0; <a class="code" href="classdgematrix.html#a67f90160391aa6dbecc09b8eb6e23246">i</a>&lt;<a class="code" href="classdgematrix.html#a60d26e3a499a235310e59a0aebbe0d0a" title="matrix column size">n</a>; <a class="code" href="classdgematrix.html#a67f90160391aa6dbecc09b8eb6e23246">i</a>++){ wr[<a class="code" href="classdgematrix.html#a67f90160391aa6dbecc09b8eb6e23246">i</a>]/=BETA[<a class="code" href="classdgematrix.html#a67f90160391aa6dbecc09b8eb6e23246">i</a>];  wi[<a class="code" href="classdgematrix.html#a67f90160391aa6dbecc09b8eb6e23246">i</a>]/=BETA[<a class="code" href="classdgematrix.html#a67f90160391aa6dbecc09b8eb6e23246">i</a>]; }
<a name="l00558"></a>00558   <span class="keyword">delete</span> [] BETA; 
<a name="l00559"></a>00559   
<a name="l00560"></a>00560   <span class="keywordflow">if</span>(INFO!=0){
<a name="l00561"></a>00561     WARNING_REPORT;
<a name="l00562"></a>00562     std::cerr &lt;&lt; <span class="stringliteral">&quot;Serious trouble happend. INFO = &quot;</span> &lt;&lt; INFO &lt;&lt; <span class="stringliteral">&quot;.&quot;</span> &lt;&lt; std::endl;
<a name="l00563"></a>00563   }
<a name="l00564"></a>00564   <span class="keywordflow">return</span> INFO;
<a name="l00565"></a>00565 }
<a name="l00566"></a>00566 
<a name="l00567"></a>00567 <span class="comment">//=============================================================================</span><span class="comment"></span>
<a name="l00568"></a>00568 <span class="comment">/*! calculate generalized eigenvalues and generalized right eigenvectors\n</span>
<a name="l00569"></a>00569 <span class="comment">  All of the arguments don&#39;t need to be initialized.</span>
<a name="l00570"></a>00570 <span class="comment">  wr, wi, vrr and vri are overwitten and become </span>
<a name="l00571"></a>00571 <span class="comment">  real and imaginary part of generalized eigenvalue </span>
<a name="l00572"></a>00572 <span class="comment">  and generalized right eigenvector, respectively. </span>
<a name="l00573"></a>00573 <span class="comment">  This matrix and matB are also overwritten.</span>
<a name="l00574"></a>00574 <span class="comment">*/</span>
<a name="l00575"></a><a class="code" href="classdgematrix.html#ae6828a6230344c9ce9463b8f229f684d">00575</a> <span class="keyword">inline</span> <span class="keywordtype">long</span> <a class="code" href="classdgematrix.html#a3374d63da801c823055f52dc7f2f55d9">dgematrix::dggev</a>(<a class="code" href="classdgematrix.html" title="Real Double-precision General Dence Matrix Class.">dgematrix</a>&amp; matB, 
<a name="l00576"></a>00576                              std::vector&lt;double&gt;&amp; wr, std::vector&lt;double&gt;&amp; wi, 
<a name="l00577"></a>00577                              std::vector&lt;dcovector&gt;&amp; vrr, 
<a name="l00578"></a>00578                              std::vector&lt;dcovector&gt;&amp; vri)
<a name="l00579"></a>00579 {VERBOSE_REPORT;
<a name="l00580"></a>00580 <span class="preprocessor">#ifdef  CPPL_DEBUG</span>
<a name="l00581"></a>00581 <span class="preprocessor"></span>  <span class="keywordflow">if</span>(<a class="code" href="classdgematrix.html#a3852aa1d0b876a74c1070e982cf388a2" title="matrix row size">m</a>!=<a class="code" href="classdgematrix.html#a60d26e3a499a235310e59a0aebbe0d0a" title="matrix column size">n</a>){
<a name="l00582"></a>00582     ERROR_REPORT;
<a name="l00583"></a>00583     std::cerr &lt;&lt; <span class="stringliteral">&quot;This matrix is not a square matrix.&quot;</span> &lt;&lt; std::endl
<a name="l00584"></a>00584               &lt;&lt; <span class="stringliteral">&quot;This matrix is (&quot;</span> &lt;&lt; <a class="code" href="classdgematrix.html#a3852aa1d0b876a74c1070e982cf388a2" title="matrix row size">m</a> &lt;&lt; <span class="stringliteral">&quot;x&quot;</span> &lt;&lt; <a class="code" href="classdgematrix.html#a60d26e3a499a235310e59a0aebbe0d0a" title="matrix column size">n</a> &lt;&lt; <span class="stringliteral">&quot;).&quot;</span> &lt;&lt; std::endl;
<a name="l00585"></a>00585     exit(1);
<a name="l00586"></a>00586   }
<a name="l00587"></a>00587   <span class="keywordflow">if</span>(matB.<a class="code" href="classdgematrix.html#a3852aa1d0b876a74c1070e982cf388a2" title="matrix row size">m</a>!=<a class="code" href="classdgematrix.html#a60d26e3a499a235310e59a0aebbe0d0a" title="matrix column size">n</a> || matB.<a class="code" href="classdgematrix.html#a60d26e3a499a235310e59a0aebbe0d0a" title="matrix column size">n</a>!=<a class="code" href="classdgematrix.html#a60d26e3a499a235310e59a0aebbe0d0a" title="matrix column size">n</a>){
<a name="l00588"></a>00588     ERROR_REPORT;
<a name="l00589"></a>00589     std::cerr &lt;&lt; <span class="stringliteral">&quot;The matrix B is not a square matrix having the same size as \&quot;this\&quot; matrix.&quot;</span> &lt;&lt; std::endl
<a name="l00590"></a>00590               &lt;&lt; <span class="stringliteral">&quot;The B matrix is (&quot;</span> &lt;&lt; matB.<a class="code" href="classdgematrix.html#a3852aa1d0b876a74c1070e982cf388a2" title="matrix row size">m</a> &lt;&lt; <span class="stringliteral">&quot;x&quot;</span> &lt;&lt; matB.<a class="code" href="classdgematrix.html#a60d26e3a499a235310e59a0aebbe0d0a" title="matrix column size">n</a> &lt;&lt; <span class="stringliteral">&quot;).&quot;</span> &lt;&lt; std::endl;
<a name="l00591"></a>00591     exit(1);
<a name="l00592"></a>00592   }
<a name="l00593"></a>00593 <span class="preprocessor">#endif//CPPL_DEBUG</span>
<a name="l00594"></a>00594 <span class="preprocessor"></span>  
<a name="l00595"></a>00595   wr.resize(<a class="code" href="classdgematrix.html#a60d26e3a499a235310e59a0aebbe0d0a" title="matrix column size">n</a>); wi.resize(<a class="code" href="classdgematrix.html#a60d26e3a499a235310e59a0aebbe0d0a" title="matrix column size">n</a>); vrr.resize(<a class="code" href="classdgematrix.html#a60d26e3a499a235310e59a0aebbe0d0a" title="matrix column size">n</a>); vri.resize(<a class="code" href="classdgematrix.html#a60d26e3a499a235310e59a0aebbe0d0a" title="matrix column size">n</a>);
<a name="l00596"></a>00596   <span class="keywordflow">for</span>(<span class="keywordtype">long</span> <a class="code" href="classdgematrix.html#a67f90160391aa6dbecc09b8eb6e23246">i</a>=0; <a class="code" href="classdgematrix.html#a67f90160391aa6dbecc09b8eb6e23246">i</a>&lt;<a class="code" href="classdgematrix.html#a60d26e3a499a235310e59a0aebbe0d0a" title="matrix column size">n</a>; <a class="code" href="classdgematrix.html#a67f90160391aa6dbecc09b8eb6e23246">i</a>++){ vrr[<a class="code" href="classdgematrix.html#a67f90160391aa6dbecc09b8eb6e23246">i</a>].resize(n); vri[<a class="code" href="classdgematrix.html#a67f90160391aa6dbecc09b8eb6e23246">i</a>].resize(n); }
<a name="l00597"></a>00597   <a class="code" href="classdgematrix.html" title="Real Double-precision General Dence Matrix Class.">dgematrix</a> VR(n,n);
<a name="l00598"></a>00598   <span class="keywordtype">char</span> JOBVL(<span class="charliteral">&#39;n&#39;</span>), JOBVR(<span class="charliteral">&#39;V&#39;</span>);
<a name="l00599"></a>00599   <span class="keywordtype">long</span> LDA(n), LDB(n), LDVL(1), LDVR(n), LWORK(8*n), INFO(1);
<a name="l00600"></a>00600   <span class="keywordtype">double</span> *BETA(<span class="keyword">new</span> <span class="keywordtype">double</span>[n]), *VL(NULL), *WORK(<span class="keyword">new</span> <span class="keywordtype">double</span>[LWORK]);
<a name="l00601"></a>00601   dggev_(JOBVL, JOBVR, n, <a class="code" href="classdgematrix.html#a006ce8c2f9add380bcb2816fdb13aa01" title="1D array to store matrix data">array</a>, LDA, matB.<a class="code" href="classdgematrix.html#a006ce8c2f9add380bcb2816fdb13aa01" title="1D array to store matrix data">array</a>, LDB, &amp;wr[0], &amp;wi[0], 
<a name="l00602"></a>00602          BETA, VL, LDVL, VR.<a class="code" href="classdgematrix.html#a006ce8c2f9add380bcb2816fdb13aa01" title="1D array to store matrix data">array</a>, LDVR, WORK, LWORK, INFO);
<a name="l00603"></a>00603   <span class="keyword">delete</span> [] WORK; <span class="keyword">delete</span> [] VL;
<a name="l00604"></a>00604   <span class="comment"></span>
<a name="l00605"></a>00605 <span class="comment">  //// reforming ////</span>
<a name="l00606"></a>00606 <span class="comment"></span>  <span class="keywordflow">for</span>(<span class="keywordtype">long</span> <a class="code" href="classdgematrix.html#a67f90160391aa6dbecc09b8eb6e23246">i</a>=0; <a class="code" href="classdgematrix.html#a67f90160391aa6dbecc09b8eb6e23246">i</a>&lt;<a class="code" href="classdgematrix.html#a60d26e3a499a235310e59a0aebbe0d0a" title="matrix column size">n</a>; <a class="code" href="classdgematrix.html#a67f90160391aa6dbecc09b8eb6e23246">i</a>++){ wr[<a class="code" href="classdgematrix.html#a67f90160391aa6dbecc09b8eb6e23246">i</a>]/=BETA[<a class="code" href="classdgematrix.html#a67f90160391aa6dbecc09b8eb6e23246">i</a>];  wi[<a class="code" href="classdgematrix.html#a67f90160391aa6dbecc09b8eb6e23246">i</a>]/=BETA[<a class="code" href="classdgematrix.html#a67f90160391aa6dbecc09b8eb6e23246">i</a>]; }
<a name="l00607"></a>00607   <span class="keyword">delete</span> [] BETA; 
<a name="l00608"></a>00608   <span class="comment"></span>
<a name="l00609"></a>00609 <span class="comment">  //// forming ////</span>
<a name="l00610"></a>00610 <span class="comment"></span>  <span class="keywordflow">for</span>(<span class="keywordtype">long</span> j=0; j&lt;<a class="code" href="classdgematrix.html#a60d26e3a499a235310e59a0aebbe0d0a" title="matrix column size">n</a>; j++){
<a name="l00611"></a>00611     <span class="keywordflow">if</span>(fabs(wi[j])&lt;DBL_MIN){
<a name="l00612"></a>00612       <span class="keywordflow">for</span>(<span class="keywordtype">long</span> <a class="code" href="classdgematrix.html#a67f90160391aa6dbecc09b8eb6e23246">i</a>=0; <a class="code" href="classdgematrix.html#a67f90160391aa6dbecc09b8eb6e23246">i</a>&lt;<a class="code" href="classdgematrix.html#a60d26e3a499a235310e59a0aebbe0d0a" title="matrix column size">n</a>; <a class="code" href="classdgematrix.html#a67f90160391aa6dbecc09b8eb6e23246">i</a>++){
<a name="l00613"></a>00613         vrr[j](<a class="code" href="classdgematrix.html#a67f90160391aa6dbecc09b8eb6e23246">i</a>) = VR(<a class="code" href="classdgematrix.html#a67f90160391aa6dbecc09b8eb6e23246">i</a>,j);
<a name="l00614"></a>00614         vri[j](<a class="code" href="classdgematrix.html#a67f90160391aa6dbecc09b8eb6e23246">i</a>) = 0.0;
<a name="l00615"></a>00615       }
<a name="l00616"></a>00616     }
<a name="l00617"></a>00617     <span class="keywordflow">else</span>{
<a name="l00618"></a>00618       <span class="keywordflow">for</span>(<span class="keywordtype">long</span> <a class="code" href="classdgematrix.html#a67f90160391aa6dbecc09b8eb6e23246">i</a>=0; <a class="code" href="classdgematrix.html#a67f90160391aa6dbecc09b8eb6e23246">i</a>&lt;<a class="code" href="classdgematrix.html#a60d26e3a499a235310e59a0aebbe0d0a" title="matrix column size">n</a>; <a class="code" href="classdgematrix.html#a67f90160391aa6dbecc09b8eb6e23246">i</a>++){
<a name="l00619"></a>00619         vrr[j](<a class="code" href="classdgematrix.html#a67f90160391aa6dbecc09b8eb6e23246">i</a>)   = VR(<a class="code" href="classdgematrix.html#a67f90160391aa6dbecc09b8eb6e23246">i</a>,j);
<a name="l00620"></a>00620         vri[j](<a class="code" href="classdgematrix.html#a67f90160391aa6dbecc09b8eb6e23246">i</a>)   = VR(<a class="code" href="classdgematrix.html#a67f90160391aa6dbecc09b8eb6e23246">i</a>,j+1);
<a name="l00621"></a>00621         vrr[j+1](<a class="code" href="classdgematrix.html#a67f90160391aa6dbecc09b8eb6e23246">i</a>) = VR(<a class="code" href="classdgematrix.html#a67f90160391aa6dbecc09b8eb6e23246">i</a>,j);
<a name="l00622"></a>00622         vri[j+1](<a class="code" href="classdgematrix.html#a67f90160391aa6dbecc09b8eb6e23246">i</a>) =-VR(<a class="code" href="classdgematrix.html#a67f90160391aa6dbecc09b8eb6e23246">i</a>,j+1);
<a name="l00623"></a>00623       }
<a name="l00624"></a>00624       j++;
<a name="l00625"></a>00625     }
<a name="l00626"></a>00626   }
<a name="l00627"></a>00627   
<a name="l00628"></a>00628   <span class="keywordflow">if</span>(INFO!=0){
<a name="l00629"></a>00629     WARNING_REPORT;
<a name="l00630"></a>00630     std::cerr &lt;&lt; <span class="stringliteral">&quot;Serious trouble happend. INFO = &quot;</span> &lt;&lt; INFO &lt;&lt; <span class="stringliteral">&quot;.&quot;</span> &lt;&lt; std::endl;
<a name="l00631"></a>00631   }
<a name="l00632"></a>00632   <span class="keywordflow">return</span> INFO;
<a name="l00633"></a>00633 }
<a name="l00634"></a>00634 
<a name="l00635"></a>00635 <span class="comment">//=============================================================================</span><span class="comment"></span>
<a name="l00636"></a>00636 <span class="comment">/*! calculate generalized eigenvalues and generalized left eigenvectors\n</span>
<a name="l00637"></a>00637 <span class="comment">  All of the arguments don&#39;t need to be initialized.</span>
<a name="l00638"></a>00638 <span class="comment">  wr, wi, vlr and vli are overwitten and become </span>
<a name="l00639"></a>00639 <span class="comment">  real and imaginary part of generalized eigenvalue </span>
<a name="l00640"></a>00640 <span class="comment">  and generalized left eigenvector, respectively. </span>
<a name="l00641"></a>00641 <span class="comment">  This matrix and matB are also overwritten.</span>
<a name="l00642"></a>00642 <span class="comment">*/</span>
<a name="l00643"></a><a class="code" href="classdgematrix.html#ac36887f04b66cde5595b5b0a8de4daf0">00643</a> <span class="keyword">inline</span> <span class="keywordtype">long</span> <a class="code" href="classdgematrix.html#a3374d63da801c823055f52dc7f2f55d9">dgematrix::dggev</a>(<a class="code" href="classdgematrix.html" title="Real Double-precision General Dence Matrix Class.">dgematrix</a>&amp; matB, 
<a name="l00644"></a>00644                              std::vector&lt;double&gt;&amp; wr, std::vector&lt;double&gt;&amp; wi, 
<a name="l00645"></a>00645                              std::vector&lt;drovector&gt;&amp; vlr, 
<a name="l00646"></a>00646                              std::vector&lt;drovector&gt;&amp; vli)
<a name="l00647"></a>00647 {VERBOSE_REPORT;
<a name="l00648"></a>00648 <span class="preprocessor">#ifdef  CPPL_DEBUG</span>
<a name="l00649"></a>00649 <span class="preprocessor"></span>  <span class="keywordflow">if</span>(<a class="code" href="classdgematrix.html#a3852aa1d0b876a74c1070e982cf388a2" title="matrix row size">m</a>!=<a class="code" href="classdgematrix.html#a60d26e3a499a235310e59a0aebbe0d0a" title="matrix column size">n</a>){
<a name="l00650"></a>00650     ERROR_REPORT;
<a name="l00651"></a>00651     std::cerr &lt;&lt; <span class="stringliteral">&quot;This matrix is not a square matrix.&quot;</span> &lt;&lt; std::endl
<a name="l00652"></a>00652               &lt;&lt; <span class="stringliteral">&quot;This matrix is (&quot;</span> &lt;&lt; <a class="code" href="classdgematrix.html#a3852aa1d0b876a74c1070e982cf388a2" title="matrix row size">m</a> &lt;&lt; <span class="stringliteral">&quot;x&quot;</span> &lt;&lt; <a class="code" href="classdgematrix.html#a60d26e3a499a235310e59a0aebbe0d0a" title="matrix column size">n</a> &lt;&lt; <span class="stringliteral">&quot;).&quot;</span> &lt;&lt; std::endl;
<a name="l00653"></a>00653     exit(1);
<a name="l00654"></a>00654   }
<a name="l00655"></a>00655   <span class="keywordflow">if</span>(matB.<a class="code" href="classdgematrix.html#a3852aa1d0b876a74c1070e982cf388a2" title="matrix row size">m</a>!=<a class="code" href="classdgematrix.html#a60d26e3a499a235310e59a0aebbe0d0a" title="matrix column size">n</a> || matB.<a class="code" href="classdgematrix.html#a60d26e3a499a235310e59a0aebbe0d0a" title="matrix column size">n</a>!=<a class="code" href="classdgematrix.html#a60d26e3a499a235310e59a0aebbe0d0a" title="matrix column size">n</a>){
<a name="l00656"></a>00656     ERROR_REPORT;
<a name="l00657"></a>00657     std::cerr &lt;&lt; <span class="stringliteral">&quot;The matrix B is not a square matrix having the same size as \&quot;this\&quot; matrix.&quot;</span> &lt;&lt; std::endl
<a name="l00658"></a>00658               &lt;&lt; <span class="stringliteral">&quot;The B matrix is (&quot;</span> &lt;&lt; matB.<a class="code" href="classdgematrix.html#a3852aa1d0b876a74c1070e982cf388a2" title="matrix row size">m</a> &lt;&lt; <span class="stringliteral">&quot;x&quot;</span> &lt;&lt; matB.<a class="code" href="classdgematrix.html#a60d26e3a499a235310e59a0aebbe0d0a" title="matrix column size">n</a> &lt;&lt; <span class="stringliteral">&quot;).&quot;</span> &lt;&lt; std::endl;
<a name="l00659"></a>00659     exit(1);
<a name="l00660"></a>00660   }
<a name="l00661"></a>00661 <span class="preprocessor">#endif//CPPL_DEBUG</span>
<a name="l00662"></a>00662 <span class="preprocessor"></span>  
<a name="l00663"></a>00663   wr.resize(<a class="code" href="classdgematrix.html#a60d26e3a499a235310e59a0aebbe0d0a" title="matrix column size">n</a>); wi.resize(<a class="code" href="classdgematrix.html#a60d26e3a499a235310e59a0aebbe0d0a" title="matrix column size">n</a>); vlr.resize(<a class="code" href="classdgematrix.html#a60d26e3a499a235310e59a0aebbe0d0a" title="matrix column size">n</a>); vli.resize(<a class="code" href="classdgematrix.html#a60d26e3a499a235310e59a0aebbe0d0a" title="matrix column size">n</a>);
<a name="l00664"></a>00664   <span class="keywordflow">for</span>(<span class="keywordtype">long</span> <a class="code" href="classdgematrix.html#a67f90160391aa6dbecc09b8eb6e23246">i</a>=0; <a class="code" href="classdgematrix.html#a67f90160391aa6dbecc09b8eb6e23246">i</a>&lt;<a class="code" href="classdgematrix.html#a60d26e3a499a235310e59a0aebbe0d0a" title="matrix column size">n</a>; <a class="code" href="classdgematrix.html#a67f90160391aa6dbecc09b8eb6e23246">i</a>++){ vlr[<a class="code" href="classdgematrix.html#a67f90160391aa6dbecc09b8eb6e23246">i</a>].resize(n); vli[<a class="code" href="classdgematrix.html#a67f90160391aa6dbecc09b8eb6e23246">i</a>].resize(n); }
<a name="l00665"></a>00665   <a class="code" href="classdgematrix.html" title="Real Double-precision General Dence Matrix Class.">dgematrix</a> VL(n,n);
<a name="l00666"></a>00666   <span class="keywordtype">char</span> JOBVL(<span class="charliteral">&#39;V&#39;</span>), JOBVR(<span class="charliteral">&#39;n&#39;</span>);
<a name="l00667"></a>00667   <span class="keywordtype">long</span> LDA(n), LDB(n), LDVL(n), LDVR(1), LWORK(8*n), INFO(1);
<a name="l00668"></a>00668   <span class="keywordtype">double</span> *BETA(<span class="keyword">new</span> <span class="keywordtype">double</span>[n]), *VR(NULL), *WORK(<span class="keyword">new</span> <span class="keywordtype">double</span>[LWORK]);
<a name="l00669"></a>00669   dggev_(JOBVL, JOBVR, n, <a class="code" href="classdgematrix.html#a006ce8c2f9add380bcb2816fdb13aa01" title="1D array to store matrix data">array</a>, LDA, matB.<a class="code" href="classdgematrix.html#a006ce8c2f9add380bcb2816fdb13aa01" title="1D array to store matrix data">array</a>, LDB, &amp;wr[0], &amp;wi[0], 
<a name="l00670"></a>00670          BETA, VL.<a class="code" href="classdgematrix.html#a006ce8c2f9add380bcb2816fdb13aa01" title="1D array to store matrix data">array</a>, LDVL, VR, LDVR, WORK, LWORK, INFO);
<a name="l00671"></a>00671   <span class="keyword">delete</span> [] WORK; <span class="keyword">delete</span> [] VR;
<a name="l00672"></a>00672   <span class="comment"></span>
<a name="l00673"></a>00673 <span class="comment">  //// reforming ////</span>
<a name="l00674"></a>00674 <span class="comment"></span>  <span class="keywordflow">for</span>(<span class="keywordtype">long</span> <a class="code" href="classdgematrix.html#a67f90160391aa6dbecc09b8eb6e23246">i</a>=0; <a class="code" href="classdgematrix.html#a67f90160391aa6dbecc09b8eb6e23246">i</a>&lt;<a class="code" href="classdgematrix.html#a60d26e3a499a235310e59a0aebbe0d0a" title="matrix column size">n</a>; <a class="code" href="classdgematrix.html#a67f90160391aa6dbecc09b8eb6e23246">i</a>++){ wr[<a class="code" href="classdgematrix.html#a67f90160391aa6dbecc09b8eb6e23246">i</a>]/=BETA[<a class="code" href="classdgematrix.html#a67f90160391aa6dbecc09b8eb6e23246">i</a>];  wi[<a class="code" href="classdgematrix.html#a67f90160391aa6dbecc09b8eb6e23246">i</a>]/=BETA[<a class="code" href="classdgematrix.html#a67f90160391aa6dbecc09b8eb6e23246">i</a>]; }
<a name="l00675"></a>00675   <span class="keyword">delete</span> [] BETA; 
<a name="l00676"></a>00676   <span class="comment"></span>
<a name="l00677"></a>00677 <span class="comment">  //// forming ////</span>
<a name="l00678"></a>00678 <span class="comment"></span>  <span class="keywordflow">for</span>(<span class="keywordtype">long</span> j=0; j&lt;<a class="code" href="classdgematrix.html#a60d26e3a499a235310e59a0aebbe0d0a" title="matrix column size">n</a>; j++){
<a name="l00679"></a>00679     <span class="keywordflow">if</span>(fabs(wi[j])&lt;DBL_MIN){
<a name="l00680"></a>00680       <span class="keywordflow">for</span>(<span class="keywordtype">long</span> <a class="code" href="classdgematrix.html#a67f90160391aa6dbecc09b8eb6e23246">i</a>=0; <a class="code" href="classdgematrix.html#a67f90160391aa6dbecc09b8eb6e23246">i</a>&lt;<a class="code" href="classdgematrix.html#a60d26e3a499a235310e59a0aebbe0d0a" title="matrix column size">n</a>; <a class="code" href="classdgematrix.html#a67f90160391aa6dbecc09b8eb6e23246">i</a>++){
<a name="l00681"></a>00681         vlr[j](<a class="code" href="classdgematrix.html#a67f90160391aa6dbecc09b8eb6e23246">i</a>) = VL(<a class="code" href="classdgematrix.html#a67f90160391aa6dbecc09b8eb6e23246">i</a>,j);
<a name="l00682"></a>00682         vli[j](<a class="code" href="classdgematrix.html#a67f90160391aa6dbecc09b8eb6e23246">i</a>) = 0.0;
<a name="l00683"></a>00683       }
<a name="l00684"></a>00684     }
<a name="l00685"></a>00685     <span class="keywordflow">else</span>{
<a name="l00686"></a>00686       <span class="keywordflow">for</span>(<span class="keywordtype">long</span> <a class="code" href="classdgematrix.html#a67f90160391aa6dbecc09b8eb6e23246">i</a>=0; <a class="code" href="classdgematrix.html#a67f90160391aa6dbecc09b8eb6e23246">i</a>&lt;<a class="code" href="classdgematrix.html#a60d26e3a499a235310e59a0aebbe0d0a" title="matrix column size">n</a>; <a class="code" href="classdgematrix.html#a67f90160391aa6dbecc09b8eb6e23246">i</a>++){
<a name="l00687"></a>00687         vlr[j](<a class="code" href="classdgematrix.html#a67f90160391aa6dbecc09b8eb6e23246">i</a>)   = VL(<a class="code" href="classdgematrix.html#a67f90160391aa6dbecc09b8eb6e23246">i</a>,j);
<a name="l00688"></a>00688         vli[j](<a class="code" href="classdgematrix.html#a67f90160391aa6dbecc09b8eb6e23246">i</a>)   =-VL(<a class="code" href="classdgematrix.html#a67f90160391aa6dbecc09b8eb6e23246">i</a>,j+1);
<a name="l00689"></a>00689         vlr[j+1](<a class="code" href="classdgematrix.html#a67f90160391aa6dbecc09b8eb6e23246">i</a>) = VL(<a class="code" href="classdgematrix.html#a67f90160391aa6dbecc09b8eb6e23246">i</a>,j);
<a name="l00690"></a>00690         vli[j+1](<a class="code" href="classdgematrix.html#a67f90160391aa6dbecc09b8eb6e23246">i</a>) = VL(<a class="code" href="classdgematrix.html#a67f90160391aa6dbecc09b8eb6e23246">i</a>,j+1);
<a name="l00691"></a>00691       }
<a name="l00692"></a>00692       j++;
<a name="l00693"></a>00693     }
<a name="l00694"></a>00694   }
<a name="l00695"></a>00695   
<a name="l00696"></a>00696   <span class="keywordflow">if</span>(INFO!=0){
<a name="l00697"></a>00697     WARNING_REPORT;
<a name="l00698"></a>00698     std::cerr &lt;&lt; <span class="stringliteral">&quot;Serious trouble happend. INFO = &quot;</span> &lt;&lt; INFO &lt;&lt; <span class="stringliteral">&quot;.&quot;</span> &lt;&lt; std::endl;
<a name="l00699"></a>00699   }
<a name="l00700"></a>00700   <span class="keywordflow">return</span> INFO;
<a name="l00701"></a>00701 }
<a name="l00702"></a>00702 <span class="comment"></span>
<a name="l00703"></a>00703 <span class="comment">///////////////////////////////////////////////////////////////////////////////</span>
<a name="l00704"></a>00704 <span class="comment">///////////////////////////////////////////////////////////////////////////////</span>
<a name="l00705"></a>00705 <span class="comment">///////////////////////////////////////////////////////////////////////////////</span>
<a name="l00706"></a>00706 <span class="comment"></span>
<a name="l00707"></a>00707 <span class="comment">//=============================================================================</span><span class="comment"></span>
<a name="l00708"></a>00708 <span class="comment">/*! compute the singular value decomposition (SVD)\n</span>
<a name="l00709"></a>00709 <span class="comment">  The argument is dgbmatrix S.</span>
<a name="l00710"></a>00710 <span class="comment">  S doesn&#39;t need to be initialized.</span>
<a name="l00711"></a>00711 <span class="comment">  S is overwitten and become singular values.</span>
<a name="l00712"></a>00712 <span class="comment">  This matrix also overwritten.</span>
<a name="l00713"></a>00713 <span class="comment">*/</span>
<a name="l00714"></a>00714 <span class="keyword">inline</span> <span class="keywordtype">long</span> <a class="code" href="classdgematrix.html#af2da13a6844756c2203bf5143d3eeb28">dgematrix::dgesvd</a>
<a name="l00715"></a><a class="code" href="classdgematrix.html#af2da13a6844756c2203bf5143d3eeb28">00715</a> (
<a name="l00716"></a>00716  <a class="code" href="classdgbmatrix.html" title="Real Double-precision General Band Matrix Class.">dgbmatrix</a>&amp; S
<a name="l00717"></a>00717  )
<a name="l00718"></a>00718 {VERBOSE_REPORT;
<a name="l00719"></a>00719   <span class="keywordtype">char</span> JOBU(<span class="charliteral">&#39;n&#39;</span>), JOBVT(<span class="charliteral">&#39;n&#39;</span>);
<a name="l00720"></a>00720   <span class="keywordtype">long</span> LDA(m), LDU(m), LDVT(n), LWORK(std::max(3*std::min(m,n)+std::max(m,n),5*std::min(m,n))), INFO(1);
<a name="l00721"></a>00721   <span class="keywordtype">double</span> *WORK(<span class="keyword">new</span> <span class="keywordtype">double</span>[LWORK]);
<a name="l00722"></a>00722   S.<a class="code" href="classdgbmatrix.html#af1e2e686717e3797cb607c86e251535c">resize</a>(m,n,0,0);
<a name="l00723"></a>00723   
<a name="l00724"></a>00724   dgesvd_(JOBU, JOBVT, m, n, array, LDA, S.<a class="code" href="classdgbmatrix.html#a3a208131351663f0aa4a95c87851d23c" title="1D array to store matrix data">array</a>, NULL, LDU, NULL, LDVT, WORK, LWORK, INFO);
<a name="l00725"></a>00725   <span class="keyword">delete</span> [] WORK;
<a name="l00726"></a>00726   
<a name="l00727"></a>00727   <span class="keywordflow">if</span>(INFO!=0){
<a name="l00728"></a>00728     WARNING_REPORT;
<a name="l00729"></a>00729     std::cerr &lt;&lt; <span class="stringliteral">&quot;Serious trouble happend. INFO = &quot;</span> &lt;&lt; INFO &lt;&lt; <span class="stringliteral">&quot;.&quot;</span> &lt;&lt; std::endl;
<a name="l00730"></a>00730   }
<a name="l00731"></a>00731   <span class="keywordflow">return</span> INFO;
<a name="l00732"></a>00732 }
<a name="l00733"></a>00733 
<a name="l00734"></a>00734 <span class="comment">//=============================================================================</span><span class="comment"></span>
<a name="l00735"></a>00735 <span class="comment">/*! compute the singular value decomposition (SVD)\n</span>
<a name="l00736"></a>00736 <span class="comment">  The arguments are dcocector S, dgematrix U and VT.</span>
<a name="l00737"></a>00737 <span class="comment">  All of them need not to be initialized.</span>
<a name="l00738"></a>00738 <span class="comment">  S, U and VT are overwitten and become singular values, </span>
<a name="l00739"></a>00739 <span class="comment">  left singular vectors,</span>
<a name="l00740"></a>00740 <span class="comment">  and right singular vectors respectively.</span>
<a name="l00741"></a>00741 <span class="comment">  This matrix also overwritten.</span>
<a name="l00742"></a>00742 <span class="comment">*/</span>
<a name="l00743"></a><a class="code" href="classdgematrix.html#af12c9079ff9bbe229a4ebfa3458e43ed">00743</a> <span class="keyword">inline</span> <span class="keywordtype">long</span> <a class="code" href="classdgematrix.html#af2da13a6844756c2203bf5143d3eeb28">dgematrix::dgesvd</a>(<a class="code" href="classdcovector.html" title="Real Double-precision Column Vector Class.">dcovector</a>&amp; S, <a class="code" href="classdgematrix.html" title="Real Double-precision General Dence Matrix Class.">dgematrix</a>&amp; U, <a class="code" href="classdgematrix.html" title="Real Double-precision General Dence Matrix Class.">dgematrix</a>&amp; VT)
<a name="l00744"></a>00744 {VERBOSE_REPORT;
<a name="l00745"></a>00745   <span class="keywordtype">char</span> JOBU(<span class="charliteral">&#39;A&#39;</span>), JOBVT(<span class="charliteral">&#39;A&#39;</span>);
<a name="l00746"></a>00746   <span class="keywordtype">long</span> LDA(<a class="code" href="classdgematrix.html#a3852aa1d0b876a74c1070e982cf388a2" title="matrix row size">m</a>), LDU(<a class="code" href="classdgematrix.html#a3852aa1d0b876a74c1070e982cf388a2" title="matrix row size">m</a>), LDVT(<a class="code" href="classdgematrix.html#a60d26e3a499a235310e59a0aebbe0d0a" title="matrix column size">n</a>),
<a name="l00747"></a>00747     LWORK(std::max(3*std::min(<a class="code" href="classdgematrix.html#a3852aa1d0b876a74c1070e982cf388a2" title="matrix row size">m</a>,<a class="code" href="classdgematrix.html#a60d26e3a499a235310e59a0aebbe0d0a" title="matrix column size">n</a>)+std::max(<a class="code" href="classdgematrix.html#a3852aa1d0b876a74c1070e982cf388a2" title="matrix row size">m</a>,<a class="code" href="classdgematrix.html#a60d26e3a499a235310e59a0aebbe0d0a" title="matrix column size">n</a>),5*std::min(<a class="code" href="classdgematrix.html#a3852aa1d0b876a74c1070e982cf388a2" title="matrix row size">m</a>,<a class="code" href="classdgematrix.html#a60d26e3a499a235310e59a0aebbe0d0a" title="matrix column size">n</a>))), INFO(1);
<a name="l00748"></a>00748   <span class="keywordtype">double</span> *WORK(<span class="keyword">new</span> <span class="keywordtype">double</span>[LWORK]);
<a name="l00749"></a>00749   S.<a class="code" href="classdcovector.html#a969d5aaa75cde3f7b1ddcfcdefddd7af">resize</a>(std::min(<a class="code" href="classdgematrix.html#a3852aa1d0b876a74c1070e982cf388a2" title="matrix row size">m</a>,<a class="code" href="classdgematrix.html#a60d26e3a499a235310e59a0aebbe0d0a" title="matrix column size">n</a>)); U.<a class="code" href="classdgematrix.html#a1e4f028034bc8b777661785e3fea6386">resize</a>(LDU,<a class="code" href="classdgematrix.html#a3852aa1d0b876a74c1070e982cf388a2" title="matrix row size">m</a>); VT.<a class="code" href="classdgematrix.html#a1e4f028034bc8b777661785e3fea6386">resize</a>(LDVT,<a class="code" href="classdgematrix.html#a60d26e3a499a235310e59a0aebbe0d0a" title="matrix column size">n</a>);
<a name="l00750"></a>00750   
<a name="l00751"></a>00751   dgesvd_(JOBU, JOBVT, <a class="code" href="classdgematrix.html#a3852aa1d0b876a74c1070e982cf388a2" title="matrix row size">m</a>, <a class="code" href="classdgematrix.html#a60d26e3a499a235310e59a0aebbe0d0a" title="matrix column size">n</a>, <a class="code" href="classdgematrix.html#a006ce8c2f9add380bcb2816fdb13aa01" title="1D array to store matrix data">array</a>, LDA, S.<a class="code" href="classdcovector.html#a9b188bf160d8f7e24ea1d41deaceaf6d" title="1D array to store vector data">array</a>, U.<a class="code" href="classdgematrix.html#a006ce8c2f9add380bcb2816fdb13aa01" title="1D array to store matrix data">array</a>,
<a name="l00752"></a>00752           LDU, VT.<a class="code" href="classdgematrix.html#a006ce8c2f9add380bcb2816fdb13aa01" title="1D array to store matrix data">array</a>, LDVT, WORK, LWORK, INFO);
<a name="l00753"></a>00753   <span class="keyword">delete</span> [] WORK;
<a name="l00754"></a>00754   
<a name="l00755"></a>00755   <span class="keywordflow">if</span>(INFO!=0){
<a name="l00756"></a>00756     WARNING_REPORT;
<a name="l00757"></a>00757     std::cerr &lt;&lt; <span class="stringliteral">&quot;Serious trouble happend. INFO = &quot;</span> &lt;&lt; INFO &lt;&lt; <span class="stringliteral">&quot;.&quot;</span> &lt;&lt; std::endl;
<a name="l00758"></a>00758   }
<a name="l00759"></a>00759   <span class="keywordflow">return</span> INFO;
<a name="l00760"></a>00760 }
<a name="l00761"></a>00761 <span class="comment"></span>
<a name="l00762"></a>00762 <span class="comment">///////////////////////////////////////////////////////////////////////////////</span>
<a name="l00763"></a>00763 <span class="comment">///////////////////////////////////////////////////////////////////////////////</span>
<a name="l00764"></a>00764 <span class="comment">///////////////////////////////////////////////////////////////////////////////</span>
<a name="l00765"></a>00765 <span class="comment"></span>
<a name="l00766"></a>00766 <span class="comment">//=============================================================================</span><span class="comment"></span>
<a name="l00767"></a>00767 <span class="comment">/*! solve the linear equality-constrained least squares (LSE) problem\n</span>
<a name="l00768"></a>00768 <span class="comment">  Input matrix and vectors, B, c, and d, are overwitten.</span>
<a name="l00769"></a>00769 <span class="comment">  This matrix is also overwritten.</span>
<a name="l00770"></a>00770 <span class="comment">  The solution vector x is to be automatically resized.</span>
<a name="l00771"></a>00771 <span class="comment">*/</span>
<a name="l00772"></a>00772 <span class="keyword">inline</span> <span class="keywordtype">long</span> <a class="code" href="classdgematrix.html#af731210d05771c998fa9ef6ddaef1bdf">dgematrix::dgglse</a>
<a name="l00773"></a><a class="code" href="classdgematrix.html#af731210d05771c998fa9ef6ddaef1bdf">00773</a> (
<a name="l00774"></a>00774  <a class="code" href="classdgematrix.html" title="Real Double-precision General Dence Matrix Class.">dgematrix</a>&amp; B,
<a name="l00775"></a>00775  <a class="code" href="classdcovector.html" title="Real Double-precision Column Vector Class.">dcovector</a>&amp; c,
<a name="l00776"></a>00776  <a class="code" href="classdcovector.html" title="Real Double-precision Column Vector Class.">dcovector</a>&amp; d,
<a name="l00777"></a>00777  <a class="code" href="classdcovector.html" title="Real Double-precision Column Vector Class.">dcovector</a>&amp; x
<a name="l00778"></a>00778  )
<a name="l00779"></a>00779 {VERBOSE_REPORT;
<a name="l00780"></a>00780 <span class="preprocessor">#ifdef  CPPL_DEBUG</span>
<a name="l00781"></a>00781 <span class="preprocessor"></span>  <span class="keywordflow">if</span>(m!=c.<a class="code" href="classdcovector.html#a08e64e5ecb5121100787d0f5d48f0af4" title="vector size">l</a>){
<a name="l00782"></a>00782     ERROR_REPORT;
<a name="l00783"></a>00783     std::cerr &lt;&lt; <span class="stringliteral">&quot;A.m and c.l should be the same.&quot;</span> &lt;&lt; std::endl
<a name="l00784"></a>00784               &lt;&lt; <span class="stringliteral">&quot;Your input was A.m=&quot;</span> &lt;&lt; m &lt;&lt; <span class="stringliteral">&quot; and c.l=&quot;</span> &lt;&lt; c.<a class="code" href="classdcovector.html#a08e64e5ecb5121100787d0f5d48f0af4" title="vector size">l</a> &lt;&lt; std::endl;
<a name="l00785"></a>00785     exit(1);
<a name="l00786"></a>00786   }
<a name="l00787"></a>00787   <span class="keywordflow">if</span>(B.<a class="code" href="classdgematrix.html#a3852aa1d0b876a74c1070e982cf388a2" title="matrix row size">m</a>!=d.<a class="code" href="classdcovector.html#a08e64e5ecb5121100787d0f5d48f0af4" title="vector size">l</a>){
<a name="l00788"></a>00788     ERROR_REPORT;
<a name="l00789"></a>00789     std::cerr &lt;&lt; <span class="stringliteral">&quot;B.m and d.l should be the same.&quot;</span> &lt;&lt; std::endl
<a name="l00790"></a>00790               &lt;&lt; <span class="stringliteral">&quot;Your input was B.m=&quot;</span> &lt;&lt; B.<a class="code" href="classdgematrix.html#a3852aa1d0b876a74c1070e982cf388a2" title="matrix row size">m</a> &lt;&lt; <span class="stringliteral">&quot; and d.l=&quot;</span> &lt;&lt; d.<a class="code" href="classdcovector.html#a08e64e5ecb5121100787d0f5d48f0af4" title="vector size">l</a> &lt;&lt; std::endl;
<a name="l00791"></a>00791     exit(1);
<a name="l00792"></a>00792   }
<a name="l00793"></a>00793   <span class="keywordflow">if</span>( !(B.<a class="code" href="classdgematrix.html#a3852aa1d0b876a74c1070e982cf388a2" title="matrix row size">m</a>&lt;=n) || !(n&lt;=m+B.<a class="code" href="classdgematrix.html#a3852aa1d0b876a74c1070e982cf388a2" title="matrix row size">m</a>) ){
<a name="l00794"></a>00794     ERROR_REPORT;
<a name="l00795"></a>00795     std::cerr &lt;&lt; <span class="stringliteral">&quot;B.m&lt;=A.n&lt;=A.m+B.m should be satisfied.&quot;</span> &lt;&lt; std::endl
<a name="l00796"></a>00796               &lt;&lt; <span class="stringliteral">&quot;Your input was B.m=&quot;</span> &lt;&lt; B.<a class="code" href="classdgematrix.html#a3852aa1d0b876a74c1070e982cf388a2" title="matrix row size">m</a> &lt;&lt; <span class="stringliteral">&quot;, A.n=&quot;</span> &lt;&lt; n &lt;&lt; <span class="stringliteral">&quot;, and A.m+B.m=&quot;</span> &lt;&lt; m+B.<a class="code" href="classdgematrix.html#a3852aa1d0b876a74c1070e982cf388a2" title="matrix row size">m</a> &lt;&lt; std::endl;
<a name="l00797"></a>00797     exit(1);
<a name="l00798"></a>00798   }
<a name="l00799"></a>00799 <span class="preprocessor">#endif//CPPL_DEBUG</span>
<a name="l00800"></a>00800 <span class="preprocessor"></span>  
<a name="l00801"></a>00801   <span class="keywordtype">long</span> lwork(-1), info(1);
<a name="l00802"></a>00802   <a class="code" href="classdcovector.html" title="Real Double-precision Column Vector Class.">dcovector</a> work(1);
<a name="l00803"></a>00803   x.<a class="code" href="classdcovector.html#a969d5aaa75cde3f7b1ddcfcdefddd7af">resize</a>(n);
<a name="l00804"></a>00804   <span class="comment"></span>
<a name="l00805"></a>00805 <span class="comment">  //////// workspace query ////////</span>
<a name="l00806"></a>00806 <span class="comment"></span>  dgglse_(m, n, B.<a class="code" href="classdgematrix.html#a3852aa1d0b876a74c1070e982cf388a2" title="matrix row size">m</a>, array, std::max(1l,m), B.<a class="code" href="classdgematrix.html#a006ce8c2f9add380bcb2816fdb13aa01" title="1D array to store matrix data">array</a>, std::max(1l,B.<a class="code" href="classdgematrix.html#a3852aa1d0b876a74c1070e982cf388a2" title="matrix row size">m</a>), c.<a class="code" href="classdcovector.html#a9b188bf160d8f7e24ea1d41deaceaf6d" title="1D array to store vector data">array</a>, d.<a class="code" href="classdcovector.html#a9b188bf160d8f7e24ea1d41deaceaf6d" title="1D array to store vector data">array</a>, x.<a class="code" href="classdcovector.html#a9b188bf160d8f7e24ea1d41deaceaf6d" title="1D array to store vector data">array</a>, work.<a class="code" href="classdcovector.html#a9b188bf160d8f7e24ea1d41deaceaf6d" title="1D array to store vector data">array</a>, lwork, info);
<a name="l00807"></a>00807   lwork =long(work(0));
<a name="l00808"></a>00808   work.<a class="code" href="classdcovector.html#a969d5aaa75cde3f7b1ddcfcdefddd7af">resize</a>(lwork);
<a name="l00809"></a>00809   info =1;
<a name="l00810"></a>00810   <span class="comment"></span>
<a name="l00811"></a>00811 <span class="comment">  //////// solve ////////</span>
<a name="l00812"></a>00812 <span class="comment"></span>  dgglse_(m, n, B.<a class="code" href="classdgematrix.html#a3852aa1d0b876a74c1070e982cf388a2" title="matrix row size">m</a>, array, std::max(1l,m), B.<a class="code" href="classdgematrix.html#a006ce8c2f9add380bcb2816fdb13aa01" title="1D array to store matrix data">array</a>, std::max(1l,B.<a class="code" href="classdgematrix.html#a3852aa1d0b876a74c1070e982cf388a2" title="matrix row size">m</a>), c.<a class="code" href="classdcovector.html#a9b188bf160d8f7e24ea1d41deaceaf6d" title="1D array to store vector data">array</a>, d.<a class="code" href="classdcovector.html#a9b188bf160d8f7e24ea1d41deaceaf6d" title="1D array to store vector data">array</a>, x.<a class="code" href="classdcovector.html#a9b188bf160d8f7e24ea1d41deaceaf6d" title="1D array to store vector data">array</a>, work.<a class="code" href="classdcovector.html#a9b188bf160d8f7e24ea1d41deaceaf6d" title="1D array to store vector data">array</a>, lwork,  info);
<a name="l00813"></a>00813   work.<a class="code" href="classdcovector.html#aa0a1e4ad0a058308b5c839609f05e973">clear</a>();
<a name="l00814"></a>00814   
<a name="l00815"></a>00815   <span class="keywordflow">if</span>(info!=0){
<a name="l00816"></a>00816     WARNING_REPORT;
<a name="l00817"></a>00817     std::cerr &lt;&lt; <span class="stringliteral">&quot;Serious trouble happend. info = &quot;</span> &lt;&lt; info &lt;&lt; <span class="stringliteral">&quot;.&quot;</span> &lt;&lt; std::endl;
<a name="l00818"></a>00818   }
<a name="l00819"></a>00819   <span class="keywordflow">return</span> info;
<a name="l00820"></a>00820 }
</pre></div></div><!-- contents -->
<!-- window showing the filter options -->
<div id="MSearchSelectWindow"
     onmouseover="return searchBox.OnSearchSelectShow()"
     onmouseout="return searchBox.OnSearchSelectHide()"
     onkeydown="return searchBox.OnSearchSelectKey(event)">
<a class="SelectItem" href="javascript:void(0)" onclick="searchBox.OnSelectItem(0)"><span class="SelectionMark">&#160;</span>All</a><a class="SelectItem" href="javascript:void(0)" onclick="searchBox.OnSelectItem(1)"><span class="SelectionMark">&#160;</span>Classes</a><a class="SelectItem" href="javascript:void(0)" onclick="searchBox.OnSelectItem(2)"><span class="SelectionMark">&#160;</span>Files</a><a class="SelectItem" href="javascript:void(0)" onclick="searchBox.OnSelectItem(3)"><span class="SelectionMark">&#160;</span>Functions</a><a class="SelectItem" href="javascript:void(0)" onclick="searchBox.OnSelectItem(4)"><span class="SelectionMark">&#160;</span>Variables</a><a class="SelectItem" href="javascript:void(0)" onclick="searchBox.OnSelectItem(5)"><span class="SelectionMark">&#160;</span>Friends</a></div>

<!-- iframe showing the search results (closed by default) -->
<div id="MSearchResultsWindow">
<iframe src="javascript:void(0)" frameborder="0" 
        name="MSearchResults" id="MSearchResults">
</iframe>
</div>



<hr class="footer"/><address class="footer"><small>
Generated by &#160;<a href="http://www.doxygen.org/index.html">
<img class="footer" src="doxygen.png" alt="doxygen"/>
</a> 1.7.6.1
</small></address>

</body>
</html>
