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 : #include <math.h>
21 : #include "solver.h"
22 :
23 :
24 0 : double** mgcLinearSystemD::NewMatrix (int N)
25 : {
26 0 : double** A = new double*[N];
27 0 : if ( !A )
28 0 : return 0;
29 :
30 0 : for (int row = 0; row < N; row++)
31 : {
32 0 : A[row] = new double[N];
33 0 : if ( !A[row] )
34 : {
35 0 : for (int i = 0; i < row; i++)
36 0 : delete[] A[i];
37 0 : delete[] A;
38 0 : return 0;
39 : }
40 0 : for (int col = 0; col < N; col++)
41 0 : A[row][col] = 0;
42 : }
43 0 : return A;
44 : }
45 :
46 0 : void mgcLinearSystemD::DeleteMatrix (int N, double** A)
47 : {
48 0 : for (int row = 0; row < N; row++)
49 0 : delete[] A[row];
50 0 : delete[] A;
51 0 : }
52 :
53 0 : double* mgcLinearSystemD::NewVector (int N)
54 : {
55 0 : double* B = new double[N];
56 0 : if ( !B )
57 0 : return 0;
58 :
59 0 : for (int row = 0; row < N; row++)
60 0 : B[row] = 0;
61 0 : return B;
62 : }
63 :
64 0 : int mgcLinearSystemD::Solve (int n, double** a, double* b)
65 : {
66 0 : int* indxc = new int[n];
67 0 : if ( !indxc )
68 0 : return 0;
69 0 : int* indxr = new int[n];
70 0 : if ( !indxr ) {
71 0 : delete[] indxc;
72 0 : return 0;
73 : }
74 0 : int* ipiv = new int[n];
75 0 : if ( !ipiv ) {
76 0 : delete[] indxc;
77 0 : delete[] indxr;
78 0 : return 0;
79 : }
80 :
81 : int i, j, k;
82 0 : int irow = 0;
83 0 : int icol = 0;
84 : double save;
85 :
86 0 : for (j = 0; j < n; j++)
87 0 : ipiv[j] = 0;
88 :
89 0 : for (i = 0; i < n; i++)
90 : {
91 0 : double big = 0;
92 0 : for (j = 0; j < n; j++)
93 : {
94 0 : if ( ipiv[j] != 1 )
95 : {
96 0 : for (k = 0; k < n; k++)
97 : {
98 0 : if ( ipiv[k] == 0 )
99 : {
100 0 : if ( fabs(a[j][k]) >= big )
101 : {
102 0 : big = fabs(a[j][k]);
103 0 : irow = j;
104 0 : icol = k;
105 : }
106 : }
107 0 : else if ( ipiv[k] > 1 )
108 : {
109 0 : delete[] ipiv;
110 0 : delete[] indxr;
111 0 : delete[] indxc;
112 0 : return 0;
113 : }
114 : }
115 : }
116 : }
117 0 : ipiv[icol]++;
118 :
119 0 : if ( irow != icol )
120 : {
121 0 : double* rowptr = a[irow];
122 0 : a[irow] = a[icol];
123 0 : a[icol] = rowptr;
124 :
125 0 : save = b[irow];
126 0 : b[irow] = b[icol];
127 0 : b[icol] = save;
128 : }
129 :
130 0 : indxr[i] = irow;
131 0 : indxc[i] = icol;
132 0 : if ( a[icol][icol] == 0 )
133 : {
134 0 : delete[] ipiv;
135 0 : delete[] indxr;
136 0 : delete[] indxc;
137 0 : return 0;
138 : }
139 :
140 0 : double pivinv = 1/a[icol][icol];
141 0 : a[icol][icol] = 1;
142 0 : for (k = 0; k < n; k++)
143 0 : a[icol][k] *= pivinv;
144 0 : b[icol] *= pivinv;
145 :
146 0 : for (j = 0; j < n; j++)
147 : {
148 0 : if ( j != icol )
149 : {
150 0 : save = a[j][icol];
151 0 : a[j][icol] = 0;
152 0 : for (k = 0; k < n; k++)
153 0 : a[j][k] -= a[icol][k]*save;
154 0 : b[j] -= b[icol]*save;
155 : }
156 : }
157 : }
158 :
159 0 : for (j = n-1; j >= 0; j--)
160 : {
161 0 : if ( indxr[j] != indxc[j] )
162 : {
163 0 : for (k = 0; k < n; k++)
164 : {
165 0 : save = a[k][indxr[j]];
166 0 : a[k][indxr[j]] = a[k][indxc[j]];
167 0 : a[k][indxc[j]] = save;
168 : }
169 : }
170 : }
171 :
172 0 : delete[] ipiv;
173 0 : delete[] indxr;
174 0 : delete[] indxc;
175 0 : return 1;
176 : }
177 :
178 : /* vim:set shiftwidth=4 softtabstop=4 expandtab: */
|
