xref: /aoo41x/main/basegfx/source/workbench/gauss.hxx (revision cdf0e10c)
1 /*************************************************************************
2  *
3  * DO NOT ALTER OR REMOVE COPYRIGHT NOTICES OR THIS FILE HEADER.
4  *
5  * Copyright 2000, 2010 Oracle and/or its affiliates.
6  *
7  * OpenOffice.org - a multi-platform office productivity suite
8  *
9  * This file is part of OpenOffice.org.
10  *
11  * OpenOffice.org is free software: you can redistribute it and/or modify
12  * it under the terms of the GNU Lesser General Public License version 3
13  * only, as published by the Free Software Foundation.
14  *
15  * OpenOffice.org is distributed in the hope that it will be useful,
16  * but WITHOUT ANY WARRANTY; without even the implied warranty of
17  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
18  * GNU Lesser General Public License version 3 for more details
19  * (a copy is included in the LICENSE file that accompanied this code).
20  *
21  * You should have received a copy of the GNU Lesser General Public License
22  * version 3 along with OpenOffice.org.  If not, see
23  * <http://www.openoffice.org/license.html>
24  * for a copy of the LGPLv3 License.
25  *
26  ************************************************************************/
27 
28 /** This method eliminates elements below main diagonal in the given
29     matrix by gaussian elimination.
30 
31     @param matrix
32     The matrix to operate on. Last column is the result vector (right
33     hand side of the linear equation). After successful termination,
34     the matrix is upper triangular. The matrix is expected to be in
35     row major order.
36 
37     @param rows
38     Number of rows in matrix
39 
40     @param cols
41     Number of columns in matrix
42 
43     @param minPivot
44     If the pivot element gets lesser than minPivot, this method fails,
45     otherwise, elimination succeeds and true is returned.
46 
47     @return true, if elimination succeeded.
48  */
49 template <class Matrix, typename BaseType>
50 bool eliminate( 	Matrix&			matrix,
51                     int				rows,
52                     int				cols,
53                     const BaseType& minPivot	)
54 {
55 	BaseType	temp;
56 	int			max, i, j, k;	/* *must* be signed, when looping like: j>=0 ! */
57 
58 	/* eliminate below main diagonal */
59 	for(i=0; i<cols-1; ++i)
60 	{
61 		/* find best pivot */
62 		max = i;
63 		for(j=i+1; j<rows; ++j)
64 			if( fabs(matrix[ j*cols + i ]) > fabs(matrix[ max*cols + i ]) )
65 				max = j;
66 
67 		/* check pivot value */
68 		if( fabs(matrix[ max*cols + i ]) < minPivot )
69 			return false;	/* pivot too small! */
70 
71 		/* interchange rows 'max' and 'i' */
72 		for(k=0; k<cols; ++k)
73 		{
74 			temp = matrix[ i*cols + k ];
75 			matrix[ i*cols + k ] = matrix[ max*cols + k ];
76 			matrix[ max*cols + k ] = temp;
77 		}
78 
79 		/* eliminate column */
80 		for(j=i+1; j<rows; ++j)
81 			for(k=cols-1; k>=i; --k)
82 				matrix[ j*cols + k ] -= matrix[ i*cols + k ] *
83 					matrix[ j*cols + i ] / matrix[ i*cols + i ];
84 	}
85 
86 	/* everything went well */
87 	return true;
88 }
89 
90 
91 /** Retrieve solution vector of linear system by substituting backwards.
92 
93 	This operation _relies_ on the previous successful
94     application of eliminate()!
95 
96     @param matrix
97     Matrix in upper diagonal form, as e.g. generated by eliminate()
98 
99     @param rows
100     Number of rows in matrix
101 
102     @param cols
103     Number of columns in matrix
104 
105     @param result
106     Result vector. Given matrix must have space for one column (rows entries).
107 
108     @return true, if back substitution was possible (i.e. no division
109     by zero occured).
110  */
111 template <class Matrix, class Vector, typename BaseType>
112 bool substitute(	const Matrix&	matrix,
113                     int				rows,
114                     int				cols,
115                     Vector&			result	)
116 {
117 	BaseType	temp;
118 	int			j,k;	/* *must* be signed, when looping like: j>=0 ! */
119 
120 	/* substitute backwards */
121 	for(j=rows-1; j>=0; --j)
122 	{
123 		temp = 0.0;
124 		for(k=j+1; k<cols-1; ++k)
125 			temp += matrix[ j*cols + k ] * result[k];
126 
127 		if( matrix[ j*cols + j ] == 0.0 )
128 			return false;	/* imminent division by zero! */
129 
130 		result[j] = (matrix[ j*cols + cols-1 ] - temp) / matrix[ j*cols + j ];
131 	}
132 
133 	/* everything went well */
134 	return true;
135 }
136 
137 
138 /** This method determines solution of given linear system, if any
139 
140 	This is a wrapper for eliminate and substitute, given matrix must
141 	contain right side of equation as the last column.
142 
143     @param matrix
144     The matrix to operate on. Last column is the result vector (right
145     hand side of the linear equation). After successful termination,
146     the matrix is upper triangular. The matrix is expected to be in
147     row major order.
148 
149     @param rows
150     Number of rows in matrix
151 
152     @param cols
153     Number of columns in matrix
154 
155     @param minPivot
156     If the pivot element gets lesser than minPivot, this method fails,
157     otherwise, elimination succeeds and true is returned.
158 
159     @return true, if elimination succeeded.
160  */
161 template <class Matrix, class Vector, typename BaseType>
162 bool solve(	Matrix&		matrix,
163             int			rows,
164             int			cols,
165             Vector&		result,
166             BaseType	minPivot	)
167 {
168 	if( eliminate<Matrix,BaseType>(matrix, rows, cols, minPivot) )
169 		return substitute<Matrix,Vector,BaseType>(matrix, rows, cols, result);
170 
171 	return false;
172 }
173