"Fossies" - the Fresh Open Source Software Archive

Member "gama-2.05/lib/gnu_gama/adj/homogenization.h" (10 May 2019, 7539 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 "homogenization.h" see the Fossies "Dox" file reference documentation.

    1 /*
    2     GNU Gama -- adjustment of geodetic networks
    3     Copyright (C) 2006  Ales Cepek <cepek@gnu.org>
    4 
    5     This file is part of the GNU Gama C++ 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 this library; if not, write to the Free Software
   19     Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
   20 */
   21 
   22 #ifndef GNU_Gama_Homogenization_gnu_gama_homogenization_homogenization_h
   23 #define GNU_Gama_Homogenization_gnu_gama_homogenization_homogenization_h
   24 
   25 
   26 #include <gnu_gama/adj/adj.h>
   27 #include <gnu_gama/adj/envelope.h>
   28 #include <matvec/covmat.h>
   29 #include <gnu_gama/size_to.h>
   30 #include <set>
   31 
   32 
   33 namespace GNU_gama {
   34 
   35 
   36   template <typename Float=double, typename Index=int>
   37   class Homogenization
   38   {
   39   public:
   40 
   41     Homogenization() : data(0), sm(0), ready(false)
   42     {
   43     }
   44     Homogenization(const AdjInputData* aid) : sm(0)
   45     {
   46       reset(aid);
   47     }
   48     ~Homogenization()
   49     {
   50       delete sm;
   51     }
   52 
   53     void reset(const AdjInputData* aid=0)
   54     {
   55       delete sm;
   56       pr.reset();
   57       sm    = 0;
   58       data  = aid;
   59       ready = false;
   60     }
   61 
   62     const SparseMatrix<Float, Index>* mat() { run(); return sm; }
   63     const Vec<Float>&                 rhs() { run(); return pr; }
   64 
   65 
   66   private:
   67 
   68     Homogenization(const Homogenization&);
   69     void operator=(const Homogenization&);
   70 
   71     const AdjInputData* data;
   72 
   73     typedef SparseMatrix<Float, Index> Sparse;
   74     typedef std::set<Index>            Indices;
   75 
   76     Sparse*        sm;
   77     Vec<Float>     pr;   // right hand side
   78     bool        ready;
   79 
   80 
   81     void run()
   82     {
   83       if (ready) return;
   84       if (!data) throw Exception::matvec(Exception::BadRank, "Homogenization : No input data");
   85 
   86       const BlockDiagonal<Float, Index>& cov = *data->cov();
   87 
   88       BlockDiagonal<Float, Index>* blockdiagonal = cov.replicate();
   89       blockdiagonal->cholDec();
   90       const BlockDiagonal<Float, Index>* bd = blockdiagonal;
   91 
   92       UpperBlockDiagonal<Float, Index> upper(bd);
   93       const Sparse* mata           = data->mat();
   94       Index total_scaled_nonzeroes = 0;
   95 
   96 
   97       /* homogenised right-hand side */
   98 
   99       pr = data->rhs();
  100       for (Index n, row=1; row<=pr.dim(); row++)   // forward substitution
  101         {
  102           const Float* b = upper.begin(row);
  103           const Float* e = upper.end  (row);
  104           const Float  x = pr(row) / *b++;
  105           pr(row) = x;
  106           n = row + 1;
  107           while(b != e)
  108             {
  109               pr(n++) -= *b++ * x;
  110             }
  111         }
  112 
  113 
  114       /* counting total number of nonzeros in scaled sparse matrix */
  115 
  116       std::vector<Index> block_cols(bd->blocks()+1);   // 1 based indexing
  117 
  118       for (Index row=1, block_index=1; block_index<=bd->blocks(); block_index++)
  119         {
  120           const Index  block_dim   = bd->dim  (block_index);
  121           const Index  block_width = bd->width(block_index);
  122 
  123           if (block_width == 0)    // uncorrelated observations
  124             {
  125               Index nonz = 0;
  126               for (Index i=1; i<=block_dim; i++, row++)
  127                 {
  128                   nonz += size_to<Index>(mata->end(row) - mata->begin(row));
  129                 }
  130               total_scaled_nonzeroes += nonz;
  131             }
  132           else                     // correlated observations
  133             {
  134               Indices indices;
  135               for (Index i=1; i<=block_dim; i++, row++)
  136                 {
  137                   const Index* n = mata->ibegin(row);
  138                   const Index* e = mata->iend  (row);
  139                   while (n != e)
  140                     {
  141                       indices.insert(*n++);
  142                     }
  143                 }
  144               block_cols[block_index] = size_to<Index>(indices.size());
  145               total_scaled_nonzeroes += block_dim*indices.size();
  146             }
  147         }
  148 
  149 
  150       sm = new Sparse(total_scaled_nonzeroes, mata->rows(), mata->columns());
  151 
  152 
  153       /* assembling scaled sparse matrix */
  154 
  155       std::vector<Index> perm(mata->columns()+1);  // block index permutation
  156 
  157       for (Index row=1, block_index=1; block_index<=bd->blocks(); block_index++)
  158         {
  159           const Float* block_b     = bd->begin(block_index);
  160           const Index  block_dim   = bd->dim  (block_index);
  161           const Index  block_width = bd->width(block_index);
  162 
  163           if (block_width == 0)    // uncorrelated observations
  164             for (Index i=1; i<=block_dim; i++, row++)
  165               {
  166                 sm->new_row();
  167                 const Float  d = *block_b++;
  168                 const Index* n = mata->ibegin(row);
  169                 const Float* b = mata->begin (row);
  170                 const Float* e = mata->end   (row);
  171                 while (b != e)
  172                   {
  173                     sm->add_element(*b++/d, *n++);
  174                   }
  175               }
  176           else                     // correlated observations
  177             {
  178               const Index bcols = block_cols[block_index];
  179               std::vector<Index> invp(bcols+1);     // block inverse permutaion
  180               Index invp_count = 0;
  181 
  182               Mat<Float> T(block_dim, bcols);       // matrix of block columns
  183               T.set_zero();
  184 
  185               /* copy block sparse columns to T */
  186 
  187               for (Index i=1; i<=block_dim; i++, row++)
  188                 {
  189                   const Float* b = mata->begin (row);
  190                   const Float* e = mata->end   (row);
  191                   const Index* n = mata->ibegin(row);
  192                   while (b != e)
  193                     {
  194                       const Index c = *n++;
  195                       if (perm[c] == 0)
  196                         {
  197                           perm[c] = ++invp_count;
  198                           invp[invp_count] = c;
  199                         }
  200 
  201                       T(i, perm[c]) = *b++;
  202                     }
  203                 }
  204 
  205               /* forward substitution for T */
  206 
  207               for (Index c=1; c<=bcols; c++)
  208                 for (Index n, r=row-block_dim, i=1; i<=block_dim; i++, r++)
  209                   {
  210                     const Float* b = upper.begin(r);
  211                     const Float* e = upper.end  (r);
  212                     const Float  x = T(i,c) / *b++;
  213                     T(i,c) = x;
  214                     n = i + 1;
  215                     while (b != e)
  216                       {
  217                         T(n++,c) -= *b++ * x;
  218                       }
  219                   }
  220 
  221               /* move transformed T to ouput sparse matrix  */
  222 
  223               for (Index i=1; i<=block_dim; i++)
  224                 {
  225                   sm->new_row();
  226                   for (Index j=1; j<=bcols; j++)
  227                     if (const Float element = T(i,j))
  228                       {
  229                         sm->add_element(element, invp[j]);
  230                       }
  231                 }
  232 
  233               for (Index i=1; i<=bcols; i++)      // clear permutation vector
  234                 {
  235                   perm[invp[i]] = 0;
  236                 }
  237             }
  238         }
  239 
  240       delete blockdiagonal;
  241       ready = true;
  242     }
  243 
  244   };
  245 
  246 }  // namespace GNU_gama
  247 
  248 #endif