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