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 "solver.h"
  22
  23
  24 double** mgcLinearSystemD::NewMatrix (int N)
  25 {
  26   double** A = new double*[N];
  27   if ( !A )
  28     return 0;
  29
  30   for (int row = 0; row < N; row++)
  31   {
  32     A[row] = new double[N];
  33     if ( !A[row] )
  34     {
  35       for (int i = 0; i < row; i++)
  36     delete[] A[i];
  37       delete[] A;
  38       return 0;
  39     }
  40     for (int col = 0; col < N; col++)
  41       A[row][col] = 0;
  42   }
  43   return A;
  44 }
  45
  46 void mgcLinearSystemD::DeleteMatrix (int N, double** A)
  47 {
  48   for (int row = 0; row < N; row++)
  49     delete[] A[row];
  50   delete[] A;
  51 }
  52
  53 double* mgcLinearSystemD::NewVector (int N)
  54 {
  55   double* B = new double[N];
  56   if ( !B )
  57     return 0;
  58
  59   for (int row = 0; row < N; row++)
  60     B[row] = 0;
  61   return B;
  62 }
  63
  64 int mgcLinearSystemD::Solve (int n, double** a, double* b)
  65 {
  66   int* indxc = new int[n];
  67   if ( !indxc )
  68     return 0;
  69   int* indxr = new int[n];
  70   if ( !indxr ) {
  71     delete[] indxc;
  72     return 0;
  73   }
  74   int* ipiv  = new int[n];
  75   if ( !ipiv ) {
  76     delete[] indxc;
  77     delete[] indxr;
  78     return 0;
  79   }
  80
  81   int i, j, k;
  82   int irow = 0;
  83   int icol = 0;
  84   double save;
  85
  86   for (j = 0; j < n; j++)
  87     ipiv[j] = 0;
  88
  89   for (i = 0; i < n; i++)
  90   {
  91     double big = 0;
  92     for (j = 0; j < n; j++)
  93     {
  94       if ( ipiv[j] != 1 )
  95       {
  96     for (k = 0; k < n; k++)
  97     {
  98       if ( ipiv[k] == 0 )
  99       {
 100         if ( fabs(a[j][k]) >= big )
 101         {
 102           big = fabs(a[j][k]);
 103           irow = j;
 104           icol = k;
 105         }
 106       }
 107       else if ( ipiv[k] > 1 )
 108       {
 109         delete[] ipiv;
 110         delete[] indxr;
 111         delete[] indxc;
 112         return 0;
 113       }
 114     }
 115       }
 116     }
 117     ipiv[icol]++;
 118
 119     if ( irow != icol )
 120     {
 121       double* rowptr = a[irow];
 122       a[irow] = a[icol];
 123       a[icol] = rowptr;
 124
 125       save = b[irow];
 126       b[irow] = b[icol];
 127       b[icol] = save;
 128     }
 129
 130     indxr[i] = irow;
 131     indxc[i] = icol;
 132     if ( a[icol][icol] == 0 )
 133     {
 134       delete[] ipiv;
 135       delete[] indxr;
 136       delete[] indxc;
 137       return 0;
 138     }
 139
 140     double pivinv = 1/a[icol][icol];
 141     a[icol][icol] = 1;
 142     for (k = 0; k < n; k++)
 143       a[icol][k] *= pivinv;
 144     b[icol] *= pivinv;
 145
 146     for (j = 0; j < n; j++)
 147     {
 148       if ( j != icol )
 149       {
 150     save = a[j][icol];
 151     a[j][icol] = 0;
 152     for (k = 0; k < n; k++)
 153       a[j][k] -= a[icol][k]*save;
 154     b[j] -= b[icol]*save;
 155       }
 156     }
 157   }
 158
 159   for (j = n-1; j >= 0; j--)
 160   {
 161     if ( indxr[j] != indxc[j] )
 162     {
 163       for (k = 0; k < n; k++)
 164       {
 165     save = a[k][indxr[j]];
 166     a[k][indxr[j]] = a[k][indxc[j]];
 167     a[k][indxc[j]] = save;
 168       }
 169     }
 170   }
 171
 172   delete[] ipiv;
 173   delete[] indxr;
 174   delete[] indxc;
 175   return 1;
 176 }
 177
 178 /* vim:set shiftwidth=4 softtabstop=4 expandtab: */