src/math/math_Householder.cxx

   1 // Copyright (c) 1997-1999 Matra Datavision
   2 // Copyright (c) 1999-2014 OPEN CASCADE SAS
   3 //
   4 // This file is part of Open CASCADE Technology software library.
   5 //
   6 // This library is free software; you can redistribute it and/or modify it under
   7 // the terms of the GNU Lesser General Public License version 2.1 as published
   8 // by the Free Software Foundation, with special exception defined in the file
   9 // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
  10 // distribution for complete text of the license and disclaimer of any warranty.
  11 //
  12 // Alternatively, this file may be used under the terms of Open CASCADE
  13 // commercial license or contractual agreement.
  14
  15 //#ifndef OCCT_DEBUG
  16 #define No_Standard_RangeError
  17 #define No_Standard_OutOfRange
  18 #define No_Standard_DimensionError
  19
  20 //#endif
  21
  22 #include <math_Householder.hxx>
  23 #include <math_Matrix.hxx>
  24 #include <Standard_ConstructionError.hxx>
  25 #include <Standard_DimensionError.hxx>
  26 #include <Standard_OutOfRange.hxx>
  27 #include <StdFail_NotDone.hxx>
  28
  29 // Cette classe decrit la methode de Householder qui transforme A en un
  30 // produit de matrice orthogonale par une triangulaire superieure. Les seconds
  31 // membres sont modifies dans le meme temps.
  32 // Les references sur le cote sont celles de l'algorithme explique en page
  33 // 90 du livre "Introduction a l'analyse numerique matricielle et a
  34 // l'optimisation." par P.G. CIARLET, edition MASSON. Les secondes
  35 // references sont celles du sous-programme HOUSEO d'Euclid.
  36 // A la difference du sous-programme Houseo, la premiere colonne n'est pas
  37 // traitee separement. Les tests effectues ont montre que le code effectue
  38 // specialement pour celle-ci etait plus long qu'une simple recopie. C'est
  39 // donc cette solution de recopie initiale qui a ete retenue.
  40 math_Householder::math_Householder(const math_Matrix& A, const math_Vector& B,
  41                                    const Standard_Real EPS):
  42                                    Sol(1, A.ColNumber(), 1, 1),
  43                                    Q(1, A.RowNumber(),
  44                                      1, A.ColNumber()) {
  45
  46   mylowerArow = A.LowerRow();
  47   mylowerAcol = A.LowerCol();
  48   myupperArow = A.UpperRow();
  49   myupperAcol = A.UpperCol();
  50   math_Matrix B1(1, B.Length(), 1, 1);
  51   B1.SetCol(1, B);
  52   Perform(A, B1, EPS);
  53 }
  54
  55
  56
  57 math_Householder::math_Householder(const math_Matrix& A, const math_Matrix& B,
  58                                    const Standard_Real EPS):
  59                                    Sol(1, A.ColNumber(),
  60                                        1, B.ColNumber()),
  61                                    Q(1, A.RowNumber(),
  62                                      A.LowerCol(), A.UpperCol()) {
  63
  64   mylowerArow = A.LowerRow();
  65   mylowerAcol = A.LowerCol();
  66   myupperArow = A.UpperRow();
  67   myupperAcol = A.UpperCol();
  68   Perform(A, B, EPS);
  69 }
  70
  71
  72 math_Householder::math_Householder(const math_Matrix& A, const math_Matrix& B,
  73                                    const Standard_Integer lowerArow,
  74                                    const Standard_Integer upperArow,
  75                                    const Standard_Integer lowerAcol,
  76                                    const Standard_Integer upperAcol,
  77                                    const Standard_Real EPS):
  78                                    Sol(1, upperAcol-lowerAcol+1,
  79                                        1, B.ColNumber()),
  80                                    Q(1, upperArow-lowerArow+1,
  81                                      1, upperAcol-lowerAcol+1) {
  82   mylowerArow = lowerArow;
  83   myupperArow = upperArow;
  84   mylowerAcol = lowerAcol;
  85   myupperAcol = upperAcol;
  86
  87   Perform(A, B, EPS);
  88 }
  89
  90
  91 void math_Householder::Perform(const math_Matrix& A, const math_Matrix& B,
  92                                const Standard_Real EPS) {
  93
  94   Standard_Integer i, j, k, n, l, m;
  95   Standard_Real scale, f, g, h = 0., alfaii;
  96   Standard_Real qki;
  97   Standard_Real cj;
  98   n = Q.ColNumber();
  99   l = Q.RowNumber();
 100   m = B.ColNumber();
 101   math_Matrix B2(1, l, 1, m);
 102
 103   Standard_Integer lbrow = B.LowerRow();
 104   for (i = 1; i <= l; i++) {
 105     for (j = 1; j <= n; j++) {
 106       Q(i, j) = A(i+mylowerArow-1, j+mylowerAcol-1);
 107     }
 108     for (j=1; j <= m; j++) {
 109       B2(i, j) = B(i+lbrow-1, j);
 110     }
 111   }
 112
 113   Standard_DimensionError_Raise_if(l != B.RowNumber() || n > l, " ");
 114
 115 // Traitement de chaque colonne de A:
 116     for (i = 1; i <= n; i++) {
 117       h = scale = 0.0;
 118         for (k = i; k <= l; k++) {
 119           qki = Q(k, i);
 120           h += qki*qki;                           // = ||a||*||a||     = EUAI
 121         }
 122         f = Q(i,i);                               // = a1              = AII
 123         g = f<1.e-15 ? Sqrt(h) : -Sqrt(h);
 124         if (fabs(g) <= EPS) {
 125           Done = Standard_False;
 126           return;
 127         }
 128         h -= f*g;                                 // = (v*v)/2         = C1
 129         alfaii = g-f;                             // = v               = ALFAII
 130         for (j =i+1; j <= n; j++) {
 131           scale = 0.0;
 132           for (k = i; k <= l; k++) {
 133            scale += Q(k,i)* Q(k,j);                //                   = SCAL
 134           }
 135           cj = (g*Q(i,j) - scale)/h;
 136           Q(i,j) = Q(i,j) - alfaii*cj;
 137           for(k= i+1; k <= l; k++) {
 138             Q(k,j) = Q(k, j) + cj * Q(k,i);
 139           }
 140         }
 141
 142 // Modification de B:
 143
 144         for (j = 1; j <= m; j++) {
 145           scale= Q(i,i)*B2(i,j);
 146           for (k = i+1; k <= l; k++) {
 147             scale += Q(k,i)*B2(k,j);
 148           }
 149           cj = (g*B2(i,j) - scale)/h;
 150           B2(i,j) = B2(i,j) - cj*alfaii;
 151           for (k = i+1; k <= l; k++) {
 152             B2(k,j) = B2(k,j) + cj*Q(k,i);
 153           }
 154         }
 155         Q(i,i) = g;
 156       }
 157
 158
 159   // Remontee:
 160   for (j = 1; j <= m; j++) {
 161     Sol(n,j) = B2(n,j)/Q(n,n);
 162     for (i = n -1; i >=1; i--) {
 163       scale= 0.0;
 164       for(k = i+1; k <= n; k++) {
 165         scale += Q(i,k) * Sol(k,j);
 166       }
 167       Sol(i,j) = (B2(i,j) - scale)/Q(i,i);
 168     }
 169   }
 170   Done = Standard_True;
 171 }
 172
 173
 174
 175
 176 void math_Householder::Dump(Standard_OStream& o) const {
 177
 178   o <<"math_Householder ";
 179    if (Done) {
 180      o << " Status = Done \n";
 181    }
 182    else {
 183      o << "Status = not Done \n";
 184    }
 185 }
 186