hwpfilter/source/solver.cxx

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