187 lines
4.1 KiB
C++
187 lines
4.1 KiB
C++
/* -*- Mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
|
|
/*************************************************************************
|
|
*
|
|
* DO NOT ALTER OR REMOVE COPYRIGHT NOTICES OR THIS FILE HEADER.
|
|
*
|
|
* Copyright 2000, 2010 Oracle and/or its affiliates.
|
|
*
|
|
* OpenOffice.org - a multi-platform office productivity suite
|
|
*
|
|
* This file is part of OpenOffice.org.
|
|
*
|
|
* OpenOffice.org is free software: you can redistribute it and/or modify
|
|
* it under the terms of the GNU Lesser General Public License version 3
|
|
* only, as published by the Free Software Foundation.
|
|
*
|
|
* OpenOffice.org is distributed in the hope that it will be useful,
|
|
* but WITHOUT ANY WARRANTY; without even the implied warranty of
|
|
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
|
* GNU Lesser General Public License version 3 for more details
|
|
* (a copy is included in the LICENSE file that accompanied this code).
|
|
*
|
|
* You should have received a copy of the GNU Lesser General Public License
|
|
* version 3 along with OpenOffice.org. If not, see
|
|
* <http://www.openoffice.org/license.html>
|
|
* for a copy of the LGPLv3 License.
|
|
*
|
|
************************************************************************/
|
|
|
|
#include <math.h>
|
|
#include "solver.h"
|
|
|
|
//---------------------------------------------------------------------------
|
|
double** mgcLinearSystemD::NewMatrix (int N)
|
|
{
|
|
double** A = new double*[N];
|
|
if ( !A )
|
|
return 0;
|
|
|
|
for (int row = 0; row < N; row++)
|
|
{
|
|
A[row] = new double[N];
|
|
if ( !A[row] )
|
|
{
|
|
for (int i = 0; i < row; i++)
|
|
delete[] A[i];
|
|
delete[] A;
|
|
return 0;
|
|
}
|
|
for (int col = 0; col < N; col++)
|
|
A[row][col] = 0;
|
|
}
|
|
return A;
|
|
}
|
|
//---------------------------------------------------------------------------
|
|
void mgcLinearSystemD::DeleteMatrix (int N, double** A)
|
|
{
|
|
for (int row = 0; row < N; row++)
|
|
delete[] A[row];
|
|
delete[] A;
|
|
}
|
|
//---------------------------------------------------------------------------
|
|
double* mgcLinearSystemD::NewVector (int N)
|
|
{
|
|
double* B = new double[N];
|
|
if ( !B )
|
|
return 0;
|
|
|
|
for (int row = 0; row < N; row++)
|
|
B[row] = 0;
|
|
return B;
|
|
}
|
|
//---------------------------------------------------------------------------
|
|
int mgcLinearSystemD::Solve (int n, double** a, double* b)
|
|
{
|
|
int* indxc = new int[n];
|
|
if ( !indxc )
|
|
return 0;
|
|
int* indxr = new int[n];
|
|
if ( !indxr ) {
|
|
delete[] indxc;
|
|
return 0;
|
|
}
|
|
int* ipiv = new int[n];
|
|
if ( !ipiv ) {
|
|
delete[] indxc;
|
|
delete[] indxr;
|
|
return 0;
|
|
}
|
|
|
|
int i, j, k;
|
|
int irow = 0;
|
|
int icol = 0;
|
|
double big, pivinv, save;
|
|
|
|
for (j = 0; j < n; j++)
|
|
ipiv[j] = 0;
|
|
|
|
for (i = 0; i < n; i++)
|
|
{
|
|
big = 0;
|
|
for (j = 0; j < n; j++)
|
|
{
|
|
if ( ipiv[j] != 1 )
|
|
{
|
|
for (k = 0; k < n; k++)
|
|
{
|
|
if ( ipiv[k] == 0 )
|
|
{
|
|
if ( fabs(a[j][k]) >= big )
|
|
{
|
|
big = fabs(a[j][k]);
|
|
irow = j;
|
|
icol = k;
|
|
}
|
|
}
|
|
else if ( ipiv[k] > 1 )
|
|
{
|
|
delete[] ipiv;
|
|
delete[] indxr;
|
|
delete[] indxc;
|
|
return 0;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
ipiv[icol]++;
|
|
|
|
if ( irow != icol )
|
|
{
|
|
double* rowptr = a[irow];
|
|
a[irow] = a[icol];
|
|
a[icol] = rowptr;
|
|
|
|
save = b[irow];
|
|
b[irow] = b[icol];
|
|
b[icol] = save;
|
|
}
|
|
|
|
indxr[i] = irow;
|
|
indxc[i] = icol;
|
|
if ( a[icol][icol] == 0 )
|
|
{
|
|
delete[] ipiv;
|
|
delete[] indxr;
|
|
delete[] indxc;
|
|
return 0;
|
|
}
|
|
|
|
pivinv = 1/a[icol][icol];
|
|
a[icol][icol] = 1;
|
|
for (k = 0; k < n; k++)
|
|
a[icol][k] *= pivinv;
|
|
b[icol] *= pivinv;
|
|
|
|
for (j = 0; j < n; j++)
|
|
{
|
|
if ( j != icol )
|
|
{
|
|
save = a[j][icol];
|
|
a[j][icol] = 0;
|
|
for (k = 0; k < n; k++)
|
|
a[j][k] -= a[icol][k]*save;
|
|
b[j] -= b[icol]*save;
|
|
}
|
|
}
|
|
}
|
|
|
|
for (j = n-1; j >= 0; j--)
|
|
{
|
|
if ( indxr[j] != indxc[j] )
|
|
{
|
|
for (k = 0; k < n; k++)
|
|
{
|
|
save = a[k][indxr[j]];
|
|
a[k][indxr[j]] = a[k][indxc[j]];
|
|
a[k][indxc[j]] = save;
|
|
}
|
|
}
|
|
}
|
|
|
|
delete[] ipiv;
|
|
delete[] indxr;
|
|
delete[] indxc;
|
|
return 1;
|
|
}
|
|
|
|
/* vim:set shiftwidth=4 softtabstop=4 expandtab: */
|