gauss.hxx

Go to the documentation of this file.

/* -*- 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 .
 */
 
#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;
}
 
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;
}
 
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: */