163 lines
5 KiB
C++
163 lines
5 KiB
C++
/* -*- Mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
|
|
/*
|
|
* This file is part of the LibreOffice project.
|
|
*
|
|
* This Source Code Form is subject to the terms of the Mozilla Public
|
|
* License, v. 2.0. If a copy of the MPL was not distributed with this
|
|
* file, You can obtain one at http://mozilla.org/MPL/2.0/.
|
|
*
|
|
* This file incorporates work covered by the following license notice:
|
|
*
|
|
* Licensed to the Apache Software Foundation (ASF) under one or more
|
|
* contributor license agreements. See the NOTICE file distributed
|
|
* with this work for additional information regarding copyright
|
|
* ownership. The ASF licenses this file to you under the Apache
|
|
* License, Version 2.0 (the "License"); you may not use this file
|
|
* except in compliance with the License. You may obtain a copy of
|
|
* the License at http://www.apache.org/licenses/LICENSE-2.0 .
|
|
*/
|
|
|
|
/** This method eliminates elements below main diagonal in the given
|
|
matrix by gaussian elimination.
|
|
|
|
@param matrix
|
|
The matrix to operate on. Last column is the result vector (right
|
|
hand side of the linear equation). After successful termination,
|
|
the matrix is upper triangular. The matrix is expected to be in
|
|
row major order.
|
|
|
|
@param rows
|
|
Number of rows in matrix
|
|
|
|
@param cols
|
|
Number of columns in matrix
|
|
|
|
@param minPivot
|
|
If the pivot element gets lesser than minPivot, this method fails,
|
|
otherwise, elimination succeeds and true is returned.
|
|
|
|
@return true, if elimination succeeded.
|
|
*/
|
|
|
|
#pragma once
|
|
|
|
template <class Matrix, typename BaseType>
|
|
bool eliminate( Matrix& matrix,
|
|
int rows,
|
|
int cols,
|
|
const BaseType& minPivot )
|
|
{
|
|
/* i, j, k *must* be signed, when looping like: j>=0 ! */
|
|
/* eliminate below main diagonal */
|
|
for(int i=0; i<cols-1; ++i)
|
|
{
|
|
/* find best pivot */
|
|
int max = i;
|
|
for(int j=i+1; j<rows; ++j)
|
|
if( fabs(matrix[ j*cols + i ]) > fabs(matrix[ max*cols + i ]) )
|
|
max = j;
|
|
|
|
/* check pivot value */
|
|
if( fabs(matrix[ max*cols + i ]) < minPivot )
|
|
return false; /* pivot too small! */
|
|
|
|
/* interchange rows 'max' and 'i' */
|
|
for(int k=0; k<cols; ++k)
|
|
{
|
|
std::swap(matrix[ i*cols + k ], matrix[ max*cols + k ]);
|
|
}
|
|
|
|
/* eliminate column */
|
|
for(int j=i+1; j<rows; ++j)
|
|
for(int k=cols-1; k>=i; --k)
|
|
matrix[ j*cols + k ] -= matrix[ i*cols + k ] *
|
|
matrix[ j*cols + i ] / matrix[ i*cols + i ];
|
|
}
|
|
|
|
/* everything went well */
|
|
return true;
|
|
}
|
|
|
|
/** Retrieve solution vector of linear system by substituting backwards.
|
|
|
|
This operation _relies_ on the previous successful
|
|
application of eliminate()!
|
|
|
|
@param matrix
|
|
Matrix in upper diagonal form, as e.g. generated by eliminate()
|
|
|
|
@param rows
|
|
Number of rows in matrix
|
|
|
|
@param cols
|
|
Number of columns in matrix
|
|
|
|
@param result
|
|
Result vector. Given matrix must have space for one column (rows entries).
|
|
|
|
@return true, if back substitution was possible (i.e. no division
|
|
by zero occurred).
|
|
*/
|
|
template <class Matrix, class Vector, typename BaseType>
|
|
bool substitute( const Matrix& matrix,
|
|
int rows,
|
|
int cols,
|
|
Vector& result )
|
|
{
|
|
BaseType temp;
|
|
|
|
/* j, k *must* be signed, when looping like: j>=0 ! */
|
|
/* substitute backwards */
|
|
for(int j=rows-1; j>=0; --j)
|
|
{
|
|
temp = 0.0;
|
|
for(int k=j+1; k<cols-1; ++k)
|
|
temp += matrix[ j*cols + k ] * result[k];
|
|
|
|
if( matrix[ j*cols + j ] == 0.0 )
|
|
return false; /* imminent division by zero! */
|
|
|
|
result[j] = (matrix[ j*cols + cols-1 ] - temp) / matrix[ j*cols + j ];
|
|
}
|
|
|
|
/* everything went well */
|
|
return true;
|
|
}
|
|
|
|
/** This method determines solution of given linear system, if any
|
|
|
|
This is a wrapper for eliminate and substitute, given matrix must
|
|
contain right side of equation as the last column.
|
|
|
|
@param matrix
|
|
The matrix to operate on. Last column is the result vector (right
|
|
hand side of the linear equation). After successful termination,
|
|
the matrix is upper triangular. The matrix is expected to be in
|
|
row major order.
|
|
|
|
@param rows
|
|
Number of rows in matrix
|
|
|
|
@param cols
|
|
Number of columns in matrix
|
|
|
|
@param minPivot
|
|
If the pivot element gets lesser than minPivot, this method fails,
|
|
otherwise, elimination succeeds and true is returned.
|
|
|
|
@return true, if elimination succeeded.
|
|
*/
|
|
template <class Matrix, class Vector, typename BaseType>
|
|
bool solve( Matrix& matrix,
|
|
int rows,
|
|
int cols,
|
|
Vector& result,
|
|
BaseType minPivot )
|
|
{
|
|
if( eliminate<Matrix,BaseType>(matrix, rows, cols, minPivot) )
|
|
return substitute<Matrix,Vector,BaseType>(matrix, rows, cols, result);
|
|
|
|
return false;
|
|
}
|
|
|
|
/* vim:set shiftwidth=4 softtabstop=4 expandtab: */
|