basegfx/source/workbench/gauss.hxx

   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 /** This method eliminates elements below main diagonal in the given
  21     matrix by gaussian elimination.
  22
  23     @param matrix
  24     The matrix to operate on. Last column is the result vector (right
  25     hand side of the linear equation). After successful termination,
  26     the matrix is upper triangular. The matrix is expected to be in
  27     row major order.
  28
  29     @param rows
  30     Number of rows in matrix
  31
  32     @param cols
  33     Number of columns in matrix
  34
  35     @param minPivot
  36     If the pivot element gets lesser than minPivot, this method fails,
  37     otherwise, elimination succeeds and true is returned.
  38
  39     @return true, if elimination succeeded.
  40  */
  41
  42 #ifndef INCLUDED_BASEGFX_SOURCE_WORKBENCH_GAUSS_HXX
  43 #define INCLUDED_BASEGFX_SOURCE_WORKBENCH_GAUSS_HXX
  44
  45 template <class Matrix, typename BaseType>
  46 bool eliminate(     Matrix&         matrix,
  47                     int             rows,
  48                     int             cols,
  49                     const BaseType& minPivot    )
  50 {
  51     BaseType    temp;
  52
  53     /* i, j, k *must* be signed, when looping like: j>=0 ! */
  54     /* eliminate below main diagonal */
  55     for(int i=0; i<cols-1; ++i)
  56     {
  57         /* find best pivot */
  58         int max = i;
  59         for(int j=i+1; j<rows; ++j)
  60             if( fabs(matrix[ j*cols + i ]) > fabs(matrix[ max*cols + i ]) )
  61                 max = j;
  62
  63         /* check pivot value */
  64         if( fabs(matrix[ max*cols + i ]) < minPivot )
  65             return false;   /* pivot too small! */
  66
  67         /* interchange rows 'max' and 'i' */
  68         for(int k=0; k<cols; ++k)
  69         {
  70             temp = matrix[ i*cols + k ];
  71             matrix[ i*cols + k ] = matrix[ max*cols + k ];
  72             matrix[ max*cols + k ] = temp;
  73         }
  74
  75         /* eliminate column */
  76         for(int j=i+1; j<rows; ++j)
  77             for(int k=cols-1; k>=i; --k)
  78                 matrix[ j*cols + k ] -= matrix[ i*cols + k ] *
  79                     matrix[ j*cols + i ] / matrix[ i*cols + i ];
  80     }
  81
  82     /* everything went well */
  83     return true;
  84 }
  85
  86 /** Retrieve solution vector of linear system by substituting backwards.
  87
  88     This operation _relies_ on the previous successful
  89     application of eliminate()!
  90
  91     @param matrix
  92     Matrix in upper diagonal form, as e.g. generated by eliminate()
  93
  94     @param rows
  95     Number of rows in matrix
  96
  97     @param cols
  98     Number of columns in matrix
  99
 100     @param result
 101     Result vector. Given matrix must have space for one column (rows entries).
 102
 103     @return true, if back substitution was possible (i.e. no division
 104     by zero occurred).
 105  */
 106 template <class Matrix, class Vector, typename BaseType>
 107 bool substitute(    const Matrix&   matrix,
 108                     int             rows,
 109                     int             cols,
 110                     Vector&         result  )
 111 {
 112     BaseType    temp;
 113
 114     /* j, k *must* be signed, when looping like: j>=0 ! */
 115     /* substitute backwards */
 116     for(int j=rows-1; j>=0; --j)
 117     {
 118         temp = 0.0;
 119         for(int k=j+1; k<cols-1; ++k)
 120             temp += matrix[ j*cols + k ] * result[k];
 121
 122         if( matrix[ j*cols + j ] == 0.0 )
 123             return false;   /* imminent division by zero! */
 124
 125         result[j] = (matrix[ j*cols + cols-1 ] - temp) / matrix[ j*cols + j ];
 126     }
 127
 128     /* everything went well */
 129     return true;
 130 }
 131
 132 /** This method determines solution of given linear system, if any
 133
 134     This is a wrapper for eliminate and substitute, given matrix must
 135     contain right side of equation as the last column.
 136
 137     @param matrix
 138     The matrix to operate on. Last column is the result vector (right
 139     hand side of the linear equation). After successful termination,
 140     the matrix is upper triangular. The matrix is expected to be in
 141     row major order.
 142
 143     @param rows
 144     Number of rows in matrix
 145
 146     @param cols
 147     Number of columns in matrix
 148
 149     @param minPivot
 150     If the pivot element gets lesser than minPivot, this method fails,
 151     otherwise, elimination succeeds and true is returned.
 152
 153     @return true, if elimination succeeded.
 154  */
 155 template <class Matrix, class Vector, typename BaseType>
 156 bool solve( Matrix&     matrix,
 157             int         rows,
 158             int         cols,
 159             Vector&     result,
 160             BaseType    minPivot    )
 161 {
 162     if( eliminate<Matrix,BaseType>(matrix, rows, cols, minPivot) )
 163         return substitute<Matrix,Vector,BaseType>(matrix, rows, cols, result);
 164
 165     return false;
 166 }
 167
 168 #endif // INCLUDED_BASEGFX_SOURCE_WORKBENCH_GAUSS_HXX
 169
 170 /* vim:set shiftwidth=4 softtabstop=4 expandtab: */