hwpfilter/source/solver.cpp

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