Line data    Source code

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