summaryrefslogtreecommitdiffstats
path: root/basegfx/source/workbench/gauss.hxx
diff options
context:
space:
mode:
authorDaniel Baumann <daniel.baumann@progress-linux.org>2024-04-27 16:51:28 +0000
committerDaniel Baumann <daniel.baumann@progress-linux.org>2024-04-27 16:51:28 +0000
commit940b4d1848e8c70ab7642901a68594e8016caffc (patch)
treeeb72f344ee6c3d9b80a7ecc079ea79e9fba8676d /basegfx/source/workbench/gauss.hxx
parentInitial commit. (diff)
downloadlibreoffice-940b4d1848e8c70ab7642901a68594e8016caffc.tar.xz
libreoffice-940b4d1848e8c70ab7642901a68594e8016caffc.zip
Adding upstream version 1:7.0.4.upstream/1%7.0.4upstream
Signed-off-by: Daniel Baumann <daniel.baumann@progress-linux.org>
Diffstat (limited to 'basegfx/source/workbench/gauss.hxx')
-rw-r--r--basegfx/source/workbench/gauss.hxx170
1 files changed, 170 insertions, 0 deletions
diff --git a/basegfx/source/workbench/gauss.hxx b/basegfx/source/workbench/gauss.hxx
new file mode 100644
index 000000000..fc352fe7e
--- /dev/null
+++ b/basegfx/source/workbench/gauss.hxx
@@ -0,0 +1,170 @@
+/* -*- 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.
+ */
+
+#ifndef INCLUDED_BASEGFX_SOURCE_WORKBENCH_GAUSS_HXX
+#define INCLUDED_BASEGFX_SOURCE_WORKBENCH_GAUSS_HXX
+
+template <class Matrix, typename BaseType>
+bool eliminate( Matrix& matrix,
+ int rows,
+ int cols,
+ const BaseType& minPivot )
+{
+ BaseType temp;
+
+ /* 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)
+ {
+ temp = matrix[ i*cols + k ];
+ matrix[ i*cols + k ] = matrix[ max*cols + k ];
+ matrix[ max*cols + k ] = temp;
+ }
+
+ /* 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;
+}
+
+#endif // INCLUDED_BASEGFX_SOURCE_WORKBENCH_GAUSS_HXX
+
+/* vim:set shiftwidth=4 softtabstop=4 expandtab: */