hwpfilter/source/solver.cpp

   1 /*************************************************************************
   2  *
   3  * DO NOT ALTER OR REMOVE COPYRIGHT NOTICES OR THIS FILE HEADER.
   4  *
   5  * Copyright 2000, 2010 Oracle and/or its affiliates.
   6  *
   7  * OpenOffice.org - a multi-platform office productivity suite
   8  *
   9  * This file is part of OpenOffice.org.
  10  *
  11  * OpenOffice.org is free software: you can redistribute it and/or modify
  12  * it under the terms of the GNU Lesser General Public License version 3
  13  * only, as published by the Free Software Foundation.
  14  *
  15  * OpenOffice.org is distributed in the hope that it will be useful,
  16  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  17  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  18  * GNU Lesser General Public License version 3 for more details
  19  * (a copy is included in the LICENSE file that accompanied this code).
  20  *
  21  * You should have received a copy of the GNU Lesser General Public License
  22  * version 3 along with OpenOffice.org.  If not, see
  23  * <http://www.openoffice.org/license.html>
  24  * for a copy of the LGPLv3 License.
  25  *
  26  ************************************************************************/
  27
  28 #include <math.h>
  29 #include "solver.h"
  30
  31 //---------------------------------------------------------------------------
  32 double** mgcLinearSystemD::NewMatrix (int N)
  33 {
  34   double** A = new double*[N];
  35   if ( !A )
  36     return 0;
  37
  38   for (int row = 0; row < N; row++)
  39   {
  40     A[row] = new double[N];
  41     if ( !A[row] )
  42     {
  43       for (int i = 0; i < row; i++)
  44     delete[] A[i];
  45       return 0;
  46     }
  47     for (int col = 0; col < N; col++)
  48       A[row][col] = 0;
  49   }
  50   return A;
  51 }
  52 //---------------------------------------------------------------------------
  53 void mgcLinearSystemD::DeleteMatrix (int N, double** A)
  54 {
  55   for (int row = 0; row < N; row++)
  56     delete[] A[row];
  57   delete[] A;
  58 }
  59 //---------------------------------------------------------------------------
  60 double* mgcLinearSystemD::NewVector (int N)
  61 {
  62   double* B = new double[N];
  63   if ( !B )
  64     return 0;
  65
  66   for (int row = 0; row < N; row++)
  67     B[row] = 0;
  68   return B;
  69 }
  70 //---------------------------------------------------------------------------
  71 int mgcLinearSystemD::Solve (int n, double** a, double* b)
  72 {
  73   int* indxc = new int[n];
  74   if ( !indxc )
  75     return 0;
  76   int* indxr = new int[n];
  77   if ( !indxr ) {
  78     delete[] indxc;
  79     return 0;
  80   }
  81   int* ipiv  = new int[n];
  82   if ( !ipiv ) {
  83     delete[] indxc;
  84     delete[] indxr;
  85     return 0;
  86   }
  87
  88   int i, j, k;
  89   int irow = 0;
  90   int icol = 0;
  91   double big, pivinv, save;
  92
  93   for (j = 0; j < n; j++)
  94     ipiv[j] = 0;
  95
  96   for (i = 0; i < n; i++)
  97   {
  98     big = 0;
  99     for (j = 0; j < n; j++)
 100     {
 101       if ( ipiv[j] != 1 )
 102       {
 103     for (k = 0; k < n; k++)
 104     {
 105       if ( ipiv[k] == 0 )
 106       {
 107         if ( fabs(a[j][k]) >= big )
 108         {
 109           big = fabs(a[j][k]);
 110           irow = j;
 111           icol = k;
 112         }
 113       }
 114       else if ( ipiv[k] > 1 )
 115       {
 116         delete[] ipiv;
 117         delete[] indxr;
 118         delete[] indxc;
 119         return 0;
 120       }
 121     }
 122       }
 123     }
 124     ipiv[icol]++;
 125
 126     if ( irow != icol )
 127     {
 128       double* rowptr = a[irow];
 129       a[irow] = a[icol];
 130       a[icol] = rowptr;
 131
 132       save = b[irow];
 133       b[irow] = b[icol];
 134       b[icol] = save;
 135     }
 136
 137     indxr[i] = irow;
 138     indxc[i] = icol;
 139     if ( a[icol][icol] == 0 )
 140     {
 141       delete[] ipiv;
 142       delete[] indxr;
 143       delete[] indxc;
 144       return 0;
 145     }
 146
 147     pivinv = 1/a[icol][icol];
 148     a[icol][icol] = 1;
 149     for (k = 0; k < n; k++)
 150       a[icol][k] *= pivinv;
 151     b[icol] *= pivinv;
 152
 153     for (j = 0; j < n; j++)
 154     {
 155       if ( j != icol )
 156       {
 157     save = a[j][icol];
 158     a[j][icol] = 0;
 159     for (k = 0; k < n; k++)
 160       a[j][k] -= a[icol][k]*save;
 161     b[j] -= b[icol]*save;
 162       }
 163     }
 164   }
 165
 166   for (j = n-1; j >= 0; j--)
 167   {
 168     if ( indxr[j] != indxc[j] )
 169     {
 170       for (k = 0; k < n; k++)
 171       {
 172     save = a[k][indxr[j]];
 173     a[k][indxr[j]] = a[k][indxc[j]];
 174     a[k][indxc[j]] = save;
 175       }
 176     }
 177   }
 178
 179   delete[] ipiv;
 180   delete[] indxr;
 181   delete[] indxc;
 182   return 1;
 183 }