hwpfilter/source/solver.cxx

   1 /* -*- Mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
   2 /*
   3  * This file is part of the LibreOffice project.
   4  *
   5  * This Source Code Form is subject to the terms of the Mozilla Public
   6  * License, v. 2.0. If a copy of the MPL was not distributed with this
   7  * file, You can obtain one at http://mozilla.org/MPL/2.0/.
   8  *
   9  * This file incorporates work covered by the following license notice:
  10  *
  11  *   Licensed to the Apache Software Foundation (ASF) under one or more
  12  *   contributor license agreements. See the NOTICE file distributed
  13  *   with this work for additional information regarding copyright
  14  *   ownership. The ASF licenses this file to you under the Apache
  15  *   License, Version 2.0 (the "License"); you may not use this file
  16  *   except in compliance with the License. You may obtain a copy of
  17  *   the License at http://www.apache.org/licenses/LICENSE-2.0 .
  18  */
  19
  20 #include <math.h>
  21 #include <memory>
  22 #include "solver.h"
  23
  24
  25 std::unique_ptr<std::unique_ptr<double[]>[]> mgcLinearSystemD::NewMatrix (int N)
  26 {
  27   std::unique_ptr<std::unique_ptr<double[]>[]> A(new std::unique_ptr<double[]>);
  28
  29   for (int row = 0; row < N; row++)
  30   {
  31     A[row].reset(new double[N]);
  32     for (int col = 0; col < N; col++)
  33       A[row][col] = 0;
  34   }
  35   return A;
  36 }
  37
  38 std::unique_ptr<double[]> mgcLinearSystemD::NewVector (int N)
  39 {
  40   std::unique_ptr<double[]> B(new double[N]);
  41
  42   for (int row = 0; row < N; row++)
  43     B[row] = 0;
  44   return B;
  45 }
  46
  47 bool mgcLinearSystemD::Solve (int n, std::unique_ptr<std::unique_ptr<double[]>[]> const & a, double* b)
  48 {
  49   std::unique_ptr<int[]> indxc( new int[n] );
  50   if ( !indxc )
  51     return false;
  52   std::unique_ptr<int[]> indxr( new int[n] );
  53   if ( !indxr ) {
  54     return false;
  55   }
  56   std::unique_ptr<int[]> ipiv( new int[n] );
  57   if ( !ipiv ) {
  58     return false;
  59   }
  60
  61   int i, j, k;
  62   int irow = 0;
  63   int icol = 0;
  64   double save;
  65
  66   for (j = 0; j < n; j++)
  67     ipiv[j] = 0;
  68
  69   for (i = 0; i < n; i++)
  70   {
  71     double big = 0;
  72     for (j = 0; j < n; j++)
  73     {
  74       if ( ipiv[j] != 1 )
  75       {
  76         for (k = 0; k < n; k++)
  77         {
  78           if ( ipiv[k] == 0 )
  79           {
  80             if ( fabs(a[j][k]) >= big )
  81             {
  82               big = fabs(a[j][k]);
  83               irow = j;
  84               icol = k;
  85             }
  86           }
  87           else if ( ipiv[k] > 1 )
  88           {
  89             return false;
  90           }
  91         }
  92       }
  93     }
  94     ipiv[icol]++;
  95
  96     if ( irow != icol )
  97     {
  98       std::swap(a[irow], a[icol]);
  99
 100       save = b[irow];
 101       b[irow] = b[icol];
 102       b[icol] = save;
 103     }
 104
 105     indxr[i] = irow;
 106     indxc[i] = icol;
 107     if ( a[icol][icol] == 0 )
 108     {
 109       return false;
 110     }
 111
 112     double pivinv = 1/a[icol][icol];
 113     a[icol][icol] = 1;
 114     for (k = 0; k < n; k++)
 115       a[icol][k] *= pivinv;
 116     b[icol] *= pivinv;
 117
 118     for (j = 0; j < n; j++)
 119     {
 120       if ( j != icol )
 121       {
 122     save = a[j][icol];
 123     a[j][icol] = 0;
 124     for (k = 0; k < n; k++)
 125       a[j][k] -= a[icol][k]*save;
 126     b[j] -= b[icol]*save;
 127       }
 128     }
 129   }
 130
 131   for (j = n-1; j >= 0; j--)
 132   {
 133     if ( indxr[j] != indxc[j] )
 134     {
 135       for (k = 0; k < n; k++)
 136       {
 137         save = a[k][indxr[j]];
 138         a[k][indxr[j]] = a[k][indxc[j]];
 139         a[k][indxc[j]] = save;
 140       }
 141     }
 142   }
 143
 144   return true;
 145 }
 146
 147 /* vim:set shiftwidth=4 softtabstop=4 expandtab: */