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   std::unique_ptr<int[]> indxr( new int[n] );
  51   std::unique_ptr<int[]> ipiv( new int[n] );
  52
  53   int i, j, k;
  54   int irow = 0;
  55   int icol = 0;
  56   double save;
  57
  58   for (j = 0; j < n; j++)
  59     ipiv[j] = 0;
  60
  61   for (i = 0; i < n; i++)
  62   {
  63     double big = 0;
  64     for (j = 0; j < n; j++)
  65     {
  66       if ( ipiv[j] != 1 )
  67       {
  68         for (k = 0; k < n; k++)
  69         {
  70           if ( ipiv[k] == 0 )
  71           {
  72             if ( fabs(a[j][k]) >= big )
  73             {
  74               big = fabs(a[j][k]);
  75               irow = j;
  76               icol = k;
  77             }
  78           }
  79           else if ( ipiv[k] > 1 )
  80           {
  81             return false;
  82           }
  83         }
  84       }
  85     }
  86     ipiv[icol]++;
  87
  88     if ( irow != icol )
  89     {
  90       std::swap(a[irow], a[icol]);
  91
  92       save = b[irow];
  93       b[irow] = b[icol];
  94       b[icol] = save;
  95     }
  96
  97     indxr[i] = irow;
  98     indxc[i] = icol;
  99     if ( a[icol][icol] == 0 )
 100     {
 101       return false;
 102     }
 103
 104     double pivinv = 1/a[icol][icol];
 105     a[icol][icol] = 1;
 106     for (k = 0; k < n; k++)
 107       a[icol][k] *= pivinv;
 108     b[icol] *= pivinv;
 109
 110     for (j = 0; j < n; j++)
 111     {
 112       if ( j != icol )
 113       {
 114         save = a[j][icol];
 115         a[j][icol] = 0;
 116         for (k = 0; k < n; k++)
 117             a[j][k] -= a[icol][k]*save;
 118         b[j] -= b[icol]*save;
 119       }
 120     }
 121   }
 122
 123   for (j = n-1; j >= 0; j--)
 124   {
 125     if ( indxr[j] != indxc[j] )
 126     {
 127       for (k = 0; k < n; k++)
 128       {
 129         save = a[k][indxr[j]];
 130         a[k][indxr[j]] = a[k][indxc[j]];
 131         a[k][indxc[j]] = save;
 132       }
 133     }
 134   }
 135
 136   return true;
 137 }
 138
 139 /* vim:set shiftwidth=4 softtabstop=4 expandtab: */