"Fossies" - the Fresh Open Source Software Archive

Member "gama-2.05/lib/matvec/svd.h" (10 May 2019, 20751 Bytes) of package /linux/privat/gama-2.05.tar.gz:


As a special service "Fossies" has tried to format the requested source page into HTML format using (guessed) C and C++ source code syntax highlighting (style: standard) with prefixed line numbers and code folding option. Alternatively you can here view or download the uninterpreted source code file. For more information about "svd.h" see the Fossies "Dox" file reference documentation.

    1 /*
    2   C++ Matrix/Vector templates (GNU Gama / matvec)
    3   Copyright (C) 1999, 2001, 2005, 2007, 2018  Ales Cepek <cepek@gnu.org>
    4 
    5   This file is part of the GNU Gama C++ Matrix/Vector template library.
    6 
    7   This library is free software; you can redistribute it and/or modify
    8   it under the terms of the GNU General Public License as published by
    9   the Free Software Foundation; either version 3 of the License, or
   10   (at your option) any later version.
   11 
   12   This library is distributed in the hope that it will be useful,
   13   but WITHOUT ANY WARRANTY; without even the implied warranty of
   14   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
   15   GNU General Public License for more details.
   16 
   17   You should have received a copy of the GNU General Public License
   18   along with GNU Gama.  If not, see <http://www.gnu.org/licenses/>.
   19 */
   20 
   21 #ifndef GNU_gama_gMatVec_MatSVD_h
   22 #define GNU_gama_gMatVec_MatSVD_h
   23 
   24 #include <cmath>
   25 #include <matvec/matvec.h>
   26 
   27 
   28 namespace GNU_gama {
   29 
   30   /** \brief Singular Value Decomposition  */
   31 
   32   /*  Singular Value Decomposition          A = U*W*trans(V)
   33       ############################
   34 
   35       SVD is based on the fortran SVD source from package CMLIB
   36 
   37       C * =====================================================================
   38       C * NIST Guide to Available Math Software.
   39       C * Source for module SVD from package CMLIB.
   40       C * Retrieved from ARNO on Thu Oct 29 08:04:40 1998.
   41       C * =====================================================================
   42       SUBROUTINE SVD(NM,M,N,A,W,MATU,U,MATV,V,IERR,RV1)
   43       C***BEGIN PROLOGUE  SVD
   44       C***REFER TO  EISDOC
   45       C
   46       C     This subroutine is a translation of the ALGOL procedure SVD,
   47       C     NUM. MATH. 14, 403-420(1970) by Golub and Reinsch.
   48       C     HANDBOOK FOR AUTO. COMP., VOL II-LINEAR ALGEBRA, 134-151(1971).
   49       C
   50       C     This subroutine determines the singular value decomposition
   51       C          T
   52       C     A=USV  of a REAL M by N rectangular matrix.  Householder
   53       C     bidiagonalization and a variant of the QR algorithm are used.
   54 
   55       -------------------------------------------------------------------------
   56 
   57       .....2011-04-17  (AC) suppressed g++ optimization (volatile variables)
   58 
   59       2001-02-30  (AC) Occasional problems with SVD convergence:
   60 
   61       Revisions marked with  #_LH_#  are made after the fortran subroutine QRBD
   62       published in: Solving Least Squares Problems by Charles L. Lawson and
   63       Richard J. Hanson, 1974 Prentice-Hall, Inc., Englewood Cliffs., N.J.,
   64       ISBN 0-13-822585-0, pp. 341 (QRBD : App. C, 298--300).
   65 
   66 
   67       2002-07-05  (AC) problems with SVD convergence:
   68 
   69       Three tests for convergence had to be rewritten to explicitly
   70       use a temporary variable s2:
   71 
   72       <    if ((s1 + ABS(rv1[L])) == s1) goto test_for_convergence;
   73       ---
   74       >    s2 = s1 + ABS(rv1[L]);
   75       >    if (s1 == s2) goto test_for_convergence;
   76 
   77       <    if (s1 + (ABS(W[L1])) == s1) break;
   78       ---
   79       >    s2 = s1 + ABS(W[L1]);
   80       >    if (s1 == s2) break;
   81 
   82       <    if (s1 + (ABS(f)) == s1) goto test_for_convergence;
   83       ---
   84       >    s2 = s1 + ABS(f);
   85       >    if (s1 == s2) goto test_for_convergence;
   86 
   87       ----------------------------------------------------------------------- */
   88 
   89   template <typename Float=double,
   90             typename Index=int,
   91             typename Exc=Exception::matvec>
   92   class SVD {
   93 
   94   public:
   95     SVD() : m(0), n(0), U_(), W_(), V_(),
   96             decomposed(0), W_tol(0), inv_W_(),  minx(all),
   97             list_min(0), minV(), U(0), V(0) {}
   98     SVD(const Mat<Float, Index, Exc>& A)
   99       : m(A.rows()), n(A.cols()), U_(A),
  100         W_(n), V_(n, n), decomposed(0), W_tol(0), inv_W_(n),
  101         minx(all), list_min(0), minV(), U(0), V(0)
  102     {
  103     }
  104     ~SVD()
  105     {
  106       delete[] list_min;
  107       delete[] U;
  108       delete[] V;
  109     }
  110 
  111     Float tol() { return W_tol; }
  112     Float tol(Float t) { W_tol = t; set_inv_W(); return W_tol; }
  113 
  114     SVD& decompose() { svd(); return *this; }
  115     SVD& reset(const Mat<Float, Index, Exc>& A);
  116     SVD& reset(const Mat<Float, Index, Exc>& A,
  117                const Vec<Float, Index, Exc>& w);
  118     SVD& clear();
  119 
  120     Float q_xx(Index, Index); // weight coefficient (xi,xj) of adj. unknowns
  121     Float q_bb(Index, Index); // weight coefficient (bi,bj) of adj. obs.
  122     Float q_bx(Index, Index); // weight coefficient (bi,xj)
  123 
  124     Index nullity()       { svd(); return defect; }
  125     bool  lindep(Index i) { svd(); return inv_W_(i)==0; } // linearly dep.
  126 
  127     void  min_x();               // minx = all
  128     void  min_x(Index, Index[]); // minx = subset
  129 
  130     const Mat<Float, Index, Exc>& SVD_U() { svd(); return U_; }
  131     const Vec<Float, Index, Exc>& SVD_W() { svd(); return W_; }
  132     const Mat<Float, Index, Exc>& SVD_V() { svd(); return V_; }
  133 
  134     void solve(const Vec<Float, Index, Exc>& rhs, Vec<Float, Index, Exc>& x);
  135 
  136   private:
  137     SVD(const SVD&);
  138     SVD& operator=(const SVD&);
  139 
  140     Index m, n;                // m = A.rows(), n = A.cols()
  141     Mat<Float, Index, Exc> U_;
  142     Vec<Float, Index, Exc> W_;
  143     Mat<Float, Index, Exc> V_;
  144 
  145     Index decomposed;
  146     void svd();
  147 
  148     volatile Float W_tol;
  149     Vec<Float, Index, Exc> inv_W_;
  150     void set_inv_W();
  151 
  152     enum { all, subset } minx;
  153     Index defect;
  154     Index  n_min;
  155     Index* list_min;
  156     Mat<Float, Index, Exc> minV;
  157     void min_subset_x();
  158 
  159     Float*  W;
  160     Float*  inv_W;
  161     Float** U;
  162     Float** V;
  163     void reset_UWV();
  164 
  165   };      /* class SVD */
  166 
  167 
  168   // ------ Exception::Singular Value Decompiosition member functions --------
  169 
  170 
  171   template <typename T> inline const T ABS(const T& x)
  172   {
  173     return (x >= T(0)) ? x : -x ;
  174   }
  175 
  176   template <typename Float, typename Index, typename Exc>
  177   void SVD<Float, Index, Exc>::set_inv_W()
  178   {
  179     if (W_tol == 0)
  180       {             // if not defined set W_tol to 1000*comp_epsilon
  181         volatile Float  eps, eps_1, eps_min, eps_max, sum;
  182         const Float one = 1;
  183 
  184         eps_min = Float();
  185         eps_max = eps = 1e-5;
  186         do
  187           {
  188             eps_1 = eps;
  189             eps = (eps_min + eps_max) / 2;
  190             sum = one + eps;
  191             if (sum == one)
  192               eps_min = eps;
  193             else
  194               eps_max = eps;
  195           } while (ABS(eps - eps_1)/eps > 0.1);
  196         W_tol = 1000*eps;
  197       }
  198 
  199     volatile Float vmax = Float();
  200     for (Index k = 1; k <= W_.dim(); k++) if (W[k] > vmax) vmax = W[k];
  201     const Float vmin = W_tol * vmax;
  202     defect = 0;
  203     for (Index i=1; i<=inv_W_.dim(); i++)
  204       if (ABS(W[i]) > vmin)
  205         inv_W[i] = 1/W[i];
  206       else
  207         {
  208           inv_W[i] = Float();
  209           defect++;
  210         }
  211 
  212   }      /* void SVD<Float, Index, Exc>::set_inv_W() */
  213 
  214 
  215   template <typename Float> inline Float PYTHAG( Float a, Float b )
  216   {
  217     volatile Float at, bt, ct;
  218 
  219     return
  220       (( at = ABS(a) )  > ( bt = ABS(b) ) )     ?
  221       (       ct = bt/at, at * std::sqrt( (Float)1 + ct * ct ) ) :
  222       (bt ? ( ct = at/bt, bt * std::sqrt( (Float)1 + ct * ct ) ) : (Float)0);
  223   }
  224 
  225 
  226   template <typename Float, typename Index, typename Exc>
  227   void SVD<Float, Index, Exc>::svd()
  228   {
  229     if (decomposed)
  230       return;
  231 
  232     const Float ZERO = 0;
  233     const Float ONE  = 1;
  234     const Float TWO  = 2;
  235 
  236     Index  i, i1, its, j, k, k1, L, L1;
  237     volatile Float  c, f, h, s, x, y, z ;
  238     volatile Float  s1=ZERO, g=ZERO, scale=ZERO, r, s2;
  239     volatile Float  tmp1;
  240     Index  mn;
  241 
  242     reset_UWV();
  243     Vec<Float, Index, Exc> rv1_(n);
  244     Float* rv1 = rv1_.begin() - 1;
  245 
  246     /* Householder reduction to bidiagonal form */
  247     for (i=1; i<=n; i++) {
  248       L = i+1;
  249       rv1[i] = scale*g ;
  250       g = s = scale = ZERO ;
  251       if (i <= m) {
  252         for (k=i; k<=m; k++) scale += ABS(U[k][i]);
  253         if (scale) {
  254           for (k=i; k<=m; k++) {
  255             tmp1 = (U[k][i] /= scale);
  256             s += tmp1*tmp1;
  257           }
  258           f = U[i][i];
  259           g = std::sqrt(s); if (f >= ZERO) g = -g;   // g = -SIGN(sqrt(s),f)
  260           h = f*g - s;
  261           U[i][i] = f - g;
  262           if (i != n)
  263             {
  264               for (j=L; j<=n; j++) {
  265                 s = ZERO;
  266                 for (k=i; k<=m; k++) s += U[k][i] * U[k][j];
  267                 f = s/h;
  268                 for (k=i; k<=m; k++) U[k][j] += f*U[k][i];
  269               }
  270             }
  271           for ( k = i ; k <= m ; k++ ) U[k][i] *= scale ;
  272         }
  273       }
  274       W[i] = scale*g;
  275       g = s = scale = ZERO;
  276       if (i <= m && i != n) {
  277         for (k=L; k<=n; k++) scale += ABS(U[i][k]);
  278         if (scale) {
  279           for (k=L; k<=n; k++) {
  280             tmp1 = (U[i][k] /= scale);
  281             s += tmp1*tmp1;
  282           }
  283           f = U[i][L];
  284           g = std::sqrt(s); if ( f >= ZERO ) g = -g; // g = -SIGN(sqrt(s),f)
  285           h = f*g - s;
  286           U[i][L] = f - g;
  287           for (k=L; k<=n; k++) rv1[k] = U[i][k] / h;
  288           if (i != m)
  289             {
  290               for (j=L; j<=m; j++) {
  291                 s = ZERO;
  292                 for (k=L; k<=n; k++) s += U[j][k] * U[i][k];
  293                 for ( k = L ; k <= n ; k++ ) U[j][k] += s * rv1[k];
  294               }
  295             }
  296           for (k=L; k<=n; k++) U[i][k] *= scale;
  297         }
  298       }
  299       // s1 = FMAX(s1, (FABS(W[i]+FABS(rv1[i])) );
  300       r = ABS(W[i]) + ABS(rv1[i]); if (r > s1) s1=r;
  301     }
  302 
  303     /* Accumulation of right-hand transformations */
  304     for (i=n; i>=1; i--) {
  305       if (i != n) {
  306         if (g) {
  307           for (j=L; j<=n; j++) // double division avoids possible underflow
  308             V[j][i] = (U[i][j] / U[i][L]) / g;
  309           for (j=L; j<=n; j++) {
  310             s = ZERO;
  311             for (k=L; k<=n; k++) s += U[i][k] * V[k][j];
  312             for (k=L; k<=n; k++) V[k][j] += s*V[k][i];
  313           }
  314         }
  315         for (j=L; j<=n; j++) V[i][j] = V[j][i] = ZERO;
  316       }
  317       V[i][i] = ONE;
  318       g = rv1[i];
  319       L = i;
  320     }
  321 
  322     /* Accumulation of left-hand transformations */
  323     mn = (m < n) ? m : n;
  324     for (i=mn; i>=1; i--) {    // i=IMIN(m,n);
  325       L = i+1;
  326       g = W[i];
  327       if (i != n)
  328         {
  329           for (j=L; j<=n; j++) U[i][j]=ZERO;
  330         }
  331       if (g)
  332         {
  333           if (i != mn)
  334             {
  335               for (j=L; j<=n; j++) {
  336                 s=ZERO;
  337                 for (k=L; k<=m; k++) s += U[k][i] * U[k][j];
  338                 f = (s / U[i][i]) / g;
  339                 for (k=i; k<=m; k++) U[k][j] += f * U[k][i];
  340               }
  341             }
  342           for (j=i; j<=m; j++) U[j][i] /= g;
  343         }
  344       else
  345         {
  346           for (j=i; j<=m; j++) U[j][i] = ZERO;
  347         }
  348       U[i][i] += ONE;
  349     }
  350 
  351     /* Diagonalization of the bidiagonal form */
  352 
  353     /* test for splitting */
  354     for (k=n; k>=1; k--)
  355       {
  356         k1  = k - 1;
  357         its = 0;
  358 
  359         /* test for splitting */
  360         for (;;)     // label 520 in fortran source
  361           {
  362             for (L=k; L>=1; L--)
  363               {
  364                 s2 = s1 + ABS(rv1[L]);
  365                 if (s1 == s2) goto test_for_convergence;
  366 
  367                 /* rv[1] is always zero, so there is no exit
  368                  * through the bottom of the loop */
  369                 L1 = L - 1;
  370                 s2 = s1 + ABS(W[L1]);
  371                 if (s1 == s2) break;
  372               }
  373 
  374             /* cancellation of rv1[L], if L greater then 1 */
  375             c = ZERO;
  376             s = ONE;
  377             for (i=L; i<=k; i++)
  378               {
  379                 f = s * rv1[i];
  380                 rv1[i] = c * rv1[i];
  381                 s2 = s1 + ABS(f);
  382                 if (s1 == s2) goto test_for_convergence;
  383                 g = W[i];
  384                 h = PYTHAG(f,g);
  385                 W[i] = h;
  386                 c =  g / h;
  387                 s = -f / h;
  388                 for (j=1; j<=m; j++)
  389                   {
  390                     y = U[j][L1];
  391                     z = U[j][i];
  392                     U[j][L1] =  y*c + z*s;
  393                     U[j][i]  = -y*s + z*c;
  394                   }
  395               }
  396 
  397           test_for_convergence:
  398 
  399             z = W[k];
  400             if (L == k) {
  401               /* W[k] is made nonnegative */
  402               if (z < ZERO)
  403                 {
  404                   W[k] = -z;
  405                   for (j=1; j<=n; j++) V[j][k] = -V[j][k];
  406                 }
  407               break;   /******   leaving for (;;) loop here   ******/
  408             }
  409 
  410             /* shift from bottom 2 by 2 minor */
  411 
  412             if ( its++ == 30)
  413               throw Exc(Exception::NoConvergence, "No convergence in SVD");
  414 
  415             x = W[L];
  416             y = W[k1];
  417             g = rv1[k1];
  418             h = rv1[k];
  419             // f = (Float)0.5*( ((g + z)/h)*((g - z)/y) + y/h -h/y);
  420             f = ((y-z)*(y+z)+(g-h)*(g+h))/(TWO*h*y);              // #_LH_#
  421             g = PYTHAG(f,ONE);
  422             s = (f >= ZERO) ? g : - g;         // SIGN(g,f)
  423             // f = x - (z/x)*z + (h/x)*(y/(f+s) - h);
  424             f = ((x-z)*(x+z) + h*(y/(f+s) - h))/x;                // #_LH_#
  425 
  426             /* next QR transformation */
  427             c = s = ONE;
  428             for (i1=L; i1<=k1; i1++) {
  429               i = i1 + 1;
  430               g = rv1[i];
  431               y = W[i];
  432               h = s*g;
  433               g = c*g;
  434               z = PYTHAG(f,h);
  435               rv1[i1] = z;
  436               c = f/z;
  437               s = h/z;
  438               f =  x*c + g*s;
  439               g = -x*s + g*c;
  440               h = y*s;
  441               y = y*c;
  442               for (j=1; j<=n; j++) {
  443                 x = V[j][i1];
  444                 z = V[j][i];
  445                 V[j][i1] =  x*c + z*s;
  446                 V[j][i]  = -x*s + z*c;
  447               }
  448               z = PYTHAG(f,h);
  449               W[i1] = z;
  450 
  451               /* rotation can be arbitrary if z is zero */
  452               if (z) {
  453                 c = f/z;
  454                 s = h/z;
  455               }
  456               f =  c*g + s*y;
  457               x = -s*g + c*y;
  458               for (j=1; j<=m; j++) {
  459                 y = U[j][i1];
  460                 z = U[j][i];
  461                 U[j][i1] =  y*c + z*s;
  462                 U[j][i]  = -y*s + z*c;
  463               }
  464             }
  465             rv1[L] = ZERO;
  466             rv1[k] = f;
  467             W[k] = x;
  468           }   // for (;;)
  469 
  470       }   // for k
  471 
  472     decomposed = 1;
  473     set_inv_W();
  474     if (defect > 0)
  475       {
  476         minV = V_;
  477         if (minx == subset) min_subset_x();
  478       }
  479 
  480   }      /* void SVD<Float, Index, Exc>::svd() */
  481 
  482 
  483   template <typename Float, typename Index, typename Exc>
  484   SVD<Float, Index, Exc>&
  485   SVD<Float, Index, Exc>::reset(const Mat<Float, Index, Exc>& A)
  486   {
  487     const Index R = A.rows();
  488     const Index C = A.cols();
  489     if (R != U_.rows() || C != U_.cols()) {
  490       m = R;
  491       n = C;
  492       W_.reset(C);
  493       V_.reset(C, C);
  494       inv_W_.reset(C);
  495     }
  496     U_ = A;
  497     reset_UWV();
  498 
  499     decomposed = 0;
  500 
  501     return *this;
  502   }      /* SVD& reset(const Mat<Float, Index, Exc>& A) */
  503 
  504 
  505   template <typename Float, typename Index, typename Exc>
  506   SVD<Float, Index, Exc>&
  507   SVD<Float, Index, Exc>::reset(const Mat<Float, Index, Exc>& A,
  508                                 const Vec<Float, Index, Exc>& w)
  509   {
  510     reset(A);
  511     volatile Float wi;
  512     for (Index i = 1; i <= A.rows(); i++)
  513       {
  514         wi = w(i);
  515         for (Index j = 1; j <= A.cols(); j++) U[i][j] *= wi;
  516       }
  517     return *this;
  518   }
  519 
  520 
  521   template <typename Float, typename Index, typename Exc>
  522   void SVD<Float, Index, Exc>::reset_UWV()
  523   {
  524     delete[] U;
  525     delete[] V;
  526 
  527     W = W_.begin() - 1;             // w[i] == W_(i)
  528     inv_W = inv_W_.begin() - 1;
  529     U = new Float*[m+1];
  530     U[1] = U_.begin() - 1;          // U[i][j] == U_(i,j)
  531     for (Index i=2; i<=m; i++)
  532       U[i] = U[i-1] + n;
  533     V = new Float*[n+1];
  534     V[1] = V_.begin() - 1;          // V[i][j] == V_(i,j)
  535     {   // for ...
  536       for (Index i=2; i<=n; i++)
  537         V[i] = V[i-1] + n;
  538     }   // for ...
  539 
  540   }      /* void SVD<Float, Index, Exc>::reset_UWV() */
  541 
  542 
  543   template <typename Float, typename Index, typename Exc>
  544   SVD<Float, Index, Exc>& SVD<Float, Index, Exc>::clear()
  545   {
  546     m = 0;
  547     n = 0;
  548     W_.reset();
  549     V_.reset();
  550     inv_W_.reset();
  551     U_.reset();
  552     minV.reset();
  553 
  554     delete[] U;
  555     delete[] V;
  556     U = V = 0;
  557 
  558     decomposed = 0;
  559 
  560     return *this;
  561   }      /* SVD& SVD<Float, Index, Exc>::clear() */
  562 
  563 
  564   template <typename Float, typename Index, typename Exc>
  565   Float SVD<Float, Index, Exc>::q_xx(Index i, Index j)
  566   {
  567     if (!(1 <= i && i <= n && 1 <= j && j <= n))
  568       throw Exc(Exception::BadRank, "Float SVD::q_xx(Index, Index)");
  569     // A = U*W*trans(V)
  570     // Covariance = V * inv_W * trans(inv_W) * trans(V);
  571     if (!decomposed) svd();
  572     volatile Float c = Float();
  573     for (Index k = 1; k <= n; k++)
  574       c += V[i][k] * inv_W[k] * inv_W[k] * V[j][k];
  575     return c;
  576   }
  577 
  578 
  579   template <typename Float, typename Index, typename Exc>
  580   Float SVD<Float, Index, Exc>::q_bb(Index i, Index j)
  581   {
  582     if (!(1 <= i && i <= m && 1 <= j && j <= m))
  583       throw Exc(Exception::BadRank, "Float SVD::q_bb(Index, Index)");
  584     // A = U*W*trans(V)
  585     // Covariance = U * trans(U);
  586     if (!decomposed) svd();
  587     volatile Float c = Float();
  588     for (Index k = 1; k <= n; k++)
  589       if (inv_W[k] != 0)
  590         c += U[i][k] * U[j][k];
  591     return c;
  592   }
  593 
  594 
  595   template <typename Float, typename Index, typename Exc>
  596   Float SVD<Float, Index, Exc>::q_bx(Index i, Index j)
  597   {
  598     if (!(1 <= i && i <= m && 1 <= j && j <= n))
  599       throw Exc(Exception::BadRank, "Float SVD::q_bx(Index, Index)");
  600     // A = U*W*trans(V)
  601     // Covariance = U * trans(V);
  602     if (!decomposed) svd();
  603     volatile Float c = Float();
  604     for (Index k = 1; k <= n; k++)
  605       c += U[i][k] * inv_W[k] * V[j][k];
  606     return c;
  607   }
  608 
  609 
  610   template <typename Float, typename Index, typename Exc>
  611   void SVD<Float, Index, Exc>::min_x()
  612   {
  613     if (decomposed && minx != all)
  614       {
  615         V_ = minV;
  616         delete[] V;
  617         V = new Float*[n+1];
  618         V[1] = V_.begin() - 1;        // V[i][j] == V_(i,j)
  619         for (Index i=2; i<=n; i++)
  620           V[i] = V[i-1] + n;
  621       }
  622     minx = all;
  623   }
  624 
  625 
  626   template <typename Float, typename Index, typename Exc>
  627   void SVD<Float, Index, Exc>::min_x(Index n, Index list[])
  628   {
  629     minx = subset;
  630     if (list_min != 0) delete[] list_min;
  631     n_min = n;
  632     list_min = new Index[n_min];
  633     for (Index i = 0; i < n; i++)
  634       list_min[i] = list[i];
  635 
  636     if (decomposed && defect != 0) {
  637       V_ =  minV;
  638       delete[] V;
  639       V = new Float*[n+1];
  640       V[1] = V_.begin() - 1;        // V[i][j] == V_(i,j)
  641       for (Index i=2; i<=n; i++)
  642         V[i] = V[i-1] + n;
  643       min_subset_x();
  644     }
  645   }      /* void SVD<Float, Index, Exc>::min_x(Index n, Index list[]) */
  646 
  647 
  648   template <typename Float, typename Index, typename Exc>
  649   void SVD<Float, Index, Exc>::min_subset_x()
  650   {
  651     using namespace std;
  652 
  653     if (defect == 0) return;
  654     if (defect > n_min)
  655       throw Exc(Exception::BadRegularization, "void SVD::min_subset_x()");
  656 
  657     Index im;
  658     volatile Float s;
  659     for (Index k = 1; k <= n; k++)
  660       if (inv_W[k] == 0)
  661         {
  662           volatile Float Vimk;
  663           s = Float();
  664           for (Index i = 0; i < n_min; i++) {
  665             im = list_min[i];
  666             Vimk = V[im][k];
  667             s += Vimk*Vimk;
  668           }
  669           s = std::sqrt(s);
  670           if (s == 0)
  671             throw Exc(Exception::BadRegularization,
  672                       "void SVD::min_subset_x()");
  673           { for (Index i = 1; i <= n; i++) V[i][k] /= s; }   // for ...
  674 
  675           for (Index j = 1; j <= n; j++)
  676             if (j != k)
  677               {
  678                 s = Float();
  679                 for (Index i = 0; i < n_min; i++) {
  680                   im = list_min[i];
  681                   s += V[im][j] * V[im][k];
  682                 }
  683                 { for (Index i = 1; i <= n; i++) V[i][j] -= s * V[i][k]; }
  684               }
  685         }
  686   }      /* void SVD<Float, Index, Exc>::min_subset_x() */
  687 
  688 
  689   template <typename Float, typename Index, typename Exc>
  690   void SVD<Float, Index, Exc>::solve(const Vec<Float, Index, Exc>& rhs,
  691                                      Vec<Float, Index, Exc>& x_)
  692   {
  693     svd();
  694     if (n != x_.dim()) x_.reset(n);
  695 
  696     Vec<Float, Index, Exc> t_(n);
  697     typename Vec<Float, Index, Exc>::const_iterator b = rhs.begin();
  698     typename Vec<Float, Index, Exc>::iterator t = t_.begin();
  699     volatile Float s;
  700 
  701     // t = trans(U)*b*inv(W);
  702     {   // for ...
  703       for (Index i=1; i<=n; ++i, ++t)
  704         {
  705           s = Float();
  706           b = rhs.begin();
  707           for (Index k=1; k<=m; k++, ++b) s += U[k][i] * (*b);
  708           *t = s * inv_W[i];
  709         }
  710     }   // for ...
  711 
  712     // x = V*t;
  713     typename Vec<Float, Index, Exc>::iterator x = x_.begin();
  714     for (Index i=1; i<=n; ++i, ++x)
  715       {
  716         s = Float();
  717         t = t_.begin();
  718         for (Index j=1; j<=n; ++j, ++t) s += V[i][j] * (*t);
  719         *x = s;
  720       }
  721   }      /* void solve(const VecVec& rhs, Vec& x_) */
  722 
  723 
  724 
  725 }   // namespace GNU_gama
  726 
  727 #endif