Skip to content

Commit 25e78ae

Browse files
dcwhitedcwhite
committed
Merge branch 'master' into issues_AC
2 parents 9616f0c + 89a1c39 commit 25e78ae

File tree

13 files changed

+4683
-4
lines changed

13 files changed

+4683
-4
lines changed

‎src/Core/Datatypes/CMakeLists.txt‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -34,6 +34,7 @@ SET(Core_Datatypes_SRCS
3434
Geometry.cc
3535
Material.cc
3636
Matrix.cc
37+
MatrixAlgorithms.cc
3738
MatrixTypeConversions.cc
3839
PropertyManagerExtensions.cc
3940
Scalar.cc
@@ -53,6 +54,7 @@ SET(Core_Datatypes_HEADERS
5354
HasId.h
5455
Material.h
5556
Matrix.h
57+
MatrixAlgorithms.h
5658
MatrixComparison.h
5759
MatrixFwd.h
5860
MatrixIO.h

‎src/Core/Datatypes/MatrixAlgorithms.cc‎

Lines changed: 305 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,305 @@
1+
/*
2+
For more information, please see: http://software.sci.utah.edu
3+
4+
The MIT License
5+
6+
Copyright (c) 2015 Scientific Computing and Imaging Institute,
7+
University of Utah.
8+
9+
10+
Permission is hereby granted, free of charge, to any person obtaining a
11+
copy of this software and associated documentation files (the "Software"),
12+
to deal in the Software without restriction, including without limitation
13+
the rights to use, copy, modify, merge, publish, distribute, sublicense,
14+
and/or sell copies of the Software, and to permit persons to whom the
15+
Software is furnished to do so, subject to the following conditions:
16+
17+
The above copyright notice and this permission notice shall be included
18+
in all copies or substantial portions of the Software.
19+
20+
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
21+
OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
22+
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
23+
THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
24+
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
25+
FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
26+
DEALINGS IN THE SOFTWARE.
27+
*/
28+
/// @todo Documentation Core/Datatypes/Legacy/Matrix/MatrixAlgorithms.cc
29+
//#include <Core/Util/Assert.h>
30+
#include <Core/Datatypes/Matrix.h>
31+
#include <Core/Datatypes/DenseMatrix.h>
32+
#include <Core/Datatypes/DenseColumnMatrix.h>
33+
#include <Core/Datatypes/SparseRowMatrix.h>
34+
#include <Core/Math/MiscMath.h>
35+
#include <Core/Datatypes/MatrixAlgorithms.h>
36+
//#include <Core/Datatypes/ColumnMatrixFunctions.h>
37+
38+
using namespace SCIRun::Core::Geometry;
39+
using namespace SCIRun::Core::Datatypes;
40+
41+
namespace SCIRun {
42+
#if SCIRUN4_TO_BE_ENABLED_LATER
43+
int
44+
MatrixAlgorithms::cg_solve(const Matrix<double>& matrix, const ColumnMatrix& rhs, ColumnMatrix& lhs)
45+
{
46+
double err;
47+
int niter;
48+
return cg_solve(matrix, rhs, lhs, err, niter);
49+
}
50+
51+
int
52+
MatrixAlgorithms::cg_solve(const Matrix<double>& matrix, const DenseMatrix& rhs, DenseMatrix& lhs)
53+
{
54+
double err;
55+
int niter;
56+
return cg_solve(matrix, rhs, lhs, err, niter);
57+
}
58+
59+
int
60+
MatrixAlgorithms::cg_solve(const Matrix<double>& matrix, const DenseMatrix& rhs, DenseMatrix& lhs,
61+
double &err, int &niter,
62+
double max_error, int toomany, int useLhsAsGuess)
63+
{
64+
if (rhs.ncols() != lhs.ncols()) return 0;
65+
for (index_type i=0; i<rhs.ncols(); i++)
66+
{
67+
ColumnMatrix rh(rhs.nrows()), lh(lhs.nrows());
68+
index_type j;
69+
for (j=0; j<rh.nrows(); j++)
70+
rh[j]=rhs[i][j];
71+
if (!cg_solve(matrix, rh, lh, err, niter, max_error,
72+
toomany, useLhsAsGuess)) return 0;
73+
for (j=0; j<rh.nrows(); j++)
74+
lhs[i][j]=lh[j];
75+
}
76+
return 1;
77+
}
78+
79+
int
80+
MatrixAlgorithms::cg_solve(const Matrix<double>& matrix, const ColumnMatrix& rhs, ColumnMatrix& lhs,
81+
double &err, int &niter,
82+
double max_error, int toomany, int useLhsAsGuess)
83+
{
84+
size_type size = matrix.nrows();
85+
niter=0;
86+
if (!useLhsAsGuess) lhs.zero();
87+
88+
if(toomany == 0) toomany=100*size;
89+
90+
ColumnMatrix diag(size), R(size), Z(size), P(size);
91+
92+
index_type i;
93+
for(i=0;i<size;i++) {
94+
if (Abs(matrix.get(i,i)>0.000001)) diag[i]=1./matrix.get(i,i);
95+
else diag[i]=1;
96+
}
97+
98+
matrix.mult(lhs, R);
99+
Sub(R, rhs, R);
100+
matrix.mult(R, Z);
101+
102+
double bnorm=rhs.vector_norm();
103+
err=R.vector_norm()/bnorm;
104+
105+
if(err == 0)
106+
{
107+
lhs=rhs;
108+
return 1;
109+
} else if (err>1000000) return 0;
110+
111+
double bkden=0;
112+
while(niter < toomany)
113+
{
114+
if(err < max_error)
115+
return 1;
116+
117+
niter++;
118+
119+
// Simple Preconditioning...
120+
Mult(Z, R, diag);
121+
122+
// Calculate coefficient bk and direction vectors p and pp
123+
double bknum=Dot(Z, R);
124+
125+
if(niter==1)
126+
{
127+
Copy(P, Z);
128+
}
129+
else
130+
{
131+
double bk=bknum/bkden;
132+
ScMult_Add(P, bk, P, Z);
133+
}
134+
135+
// Calculate coefficient ak, new iterate x and new residuals r and rr
136+
matrix.mult(P, Z);
137+
bkden=bknum;
138+
double akden=Dot(Z, P);
139+
140+
double ak=bknum/akden;
141+
ScMult_Add(lhs, ak, P, lhs);
142+
ScMult_Add(R, -ak, Z, R);
143+
144+
err=R.vector_norm()/bnorm;
145+
if (err>1000000) return 0;
146+
}
147+
return 0;
148+
}
149+
150+
int
151+
MatrixAlgorithms::bicg_solve(const Matrix<double>& matrix, const ColumnMatrix& rhs, ColumnMatrix& lhs)
152+
{
153+
double err;
154+
int niter;
155+
return bicg_solve(matrix, rhs, lhs, err, niter);
156+
}
157+
158+
int
159+
MatrixAlgorithms::bicg_solve(const Matrix<double>& matrix, const DenseMatrix& rhs, DenseMatrix& lhs)
160+
{
161+
double err;
162+
int niter;
163+
return bicg_solve(matrix, rhs, lhs, err, niter);
164+
}
165+
166+
int
167+
MatrixAlgorithms::bicg_solve(const Matrix<double>& matrix, const DenseMatrix& rhs, DenseMatrix& lhs,
168+
double &err, int &niter,
169+
double max_error, int /*toomany*/, int useLhsAsGuess)
170+
{
171+
if (rhs.ncols() != lhs.ncols()) return 0;
172+
for (index_type i=0; i<rhs.ncols(); i++)
173+
{
174+
ColumnMatrix rh(rhs.nrows()), lh(lhs.nrows());
175+
index_type j;
176+
for (j=0; j<rh.nrows(); j++)
177+
rh[j]=rhs[i][j];
178+
if (!bicg_solve(matrix, rh, lh, err, niter, max_error, useLhsAsGuess)) return 0;
179+
for (j=0; j<rh.nrows(); j++)
180+
lhs[i][j]=lh[j];
181+
}
182+
return 1;
183+
}
184+
185+
186+
int
187+
MatrixAlgorithms::bicg_solve(const Matrix<double>& matrix, const ColumnMatrix& rhs, ColumnMatrix& lhs,
188+
double &err, int &niter,
189+
double max_error, int toomany, int useLhsAsGuess)
190+
{
191+
size_type size=matrix.nrows();
192+
niter=0;
193+
if (!useLhsAsGuess) lhs.zero();
194+
195+
if(toomany == 0) toomany=100*size;
196+
197+
ColumnMatrix diag(size), R(size), R1(size), Z(size), Z1(size),
198+
P(size), P1(size);
199+
200+
index_type i;
201+
for(i=0;i<size;i++)
202+
{
203+
if (Abs(matrix.get(i,i)>0.000001)) diag[i]=1./matrix.get(i,i);
204+
else diag[i]=1;
205+
}
206+
207+
208+
matrix.mult(lhs, R );
209+
Sub(R, rhs, R );
210+
211+
double bnorm=rhs.vector_norm();
212+
err=R.vector_norm()/bnorm;
213+
214+
if(err == 0)
215+
{
216+
lhs=rhs;
217+
return 1;
218+
}
219+
else
220+
{
221+
if (err>1000000) return 0;
222+
}
223+
224+
// BiCG
225+
Copy(R1, R);
226+
227+
double bkden=0;
228+
while(niter < toomany)
229+
{
230+
if(err < max_error)
231+
return 1;
232+
233+
niter++;
234+
235+
// Simple Preconditioning...
236+
Mult(Z, R, diag);
237+
// BiCG
238+
Mult(Z1, R1, diag);
239+
240+
// Calculate coefficient bk and direction vectors p and pp
241+
// BiCG - change R->R1
242+
double bknum=Dot(Z, R1);
243+
244+
// BiCG
245+
if ( bknum == 0 )
246+
{
247+
return 1;
248+
}
249+
250+
if(niter==1)
251+
{
252+
Copy(P, Z);
253+
// BiCG
254+
Copy(P1, Z1);
255+
}
256+
else
257+
{
258+
double bk=bknum/bkden;
259+
ScMult_Add(P, bk, P, Z);
260+
// BiCG
261+
ScMult_Add(P1, bk, P1, Z1);
262+
}
263+
264+
// Calculate coefficient ak, new iterate x and new residuals r and rr
265+
matrix.mult(P, Z);
266+
bkden=bknum;
267+
268+
// BiCG
269+
matrix.mult_transpose(P1, Z1);
270+
271+
// BiCG = change P -> P1
272+
double akden=Dot(Z, P1);
273+
274+
double ak=bknum/akden;
275+
ScMult_Add(lhs, ak, P, lhs);
276+
ScMult_Add(R, -ak, Z, R);
277+
// BiCG
278+
ScMult_Add(R1, -ak, Z1, R1);
279+
280+
err=R.vector_norm()/bnorm;
281+
282+
if (err>1000000) return 0;
283+
}
284+
return 0;
285+
}
286+
#endif
287+
288+
Transform
289+
MatrixAlgorithms::matrix_to_transform(const Matrix& matrix)
290+
{
291+
Transform t;
292+
if (matrix.nrows() != 4 || matrix.ncols() != 4)
293+
{
294+
std::cerr << "Error: transform matrix must be 4x4.\n";
295+
return t;
296+
}
297+
double dummy[16];
298+
int cnt=0;
299+
for (index_type i=0; i<4; i++)
300+
for (index_type j=0; j<4; j++, cnt++)
301+
dummy[cnt] = matrix.get(i,j);
302+
t.set(dummy);
303+
return t;
304+
}
305+
} // End namespace SCIRun

‎src/Core/Datatypes/MatrixAlgorithms.h‎

Lines changed: 69 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,69 @@
1+
/*
2+
For more information, please see: http://software.sci.utah.edu
3+
4+
The MIT License
5+
6+
Copyright (c) 2015 Scientific Computing and Imaging Institute,
7+
University of Utah.
8+
9+
10+
Permission is hereby granted, free of charge, to any person obtaining a
11+
copy of this software and associated documentation files (the "Software"),
12+
to deal in the Software without restriction, including without limitation
13+
the rights to use, copy, modify, merge, publish, distribute, sublicense,
14+
and/or sell copies of the Software, and to permit persons to whom the
15+
Software is furnished to do so, subject to the following conditions:
16+
17+
The above copyright notice and this permission notice shall be included
18+
in all copies or substantial portions of the Software.
19+
20+
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
21+
OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
22+
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
23+
THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
24+
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
25+
FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
26+
DEALINGS IN THE SOFTWARE.
27+
*/
28+
/// @todo Documentation Core/Datatypes/Legacy/Matrix/MatrixAlgorithms.h
29+
#ifndef CORE_DATATYPES_MATRIX_ALGORITHMS_H
30+
#define CORE_DATATYPES_MATRIX_ALGORITHMS_H
31+
32+
#include <Core/Datatypes/MatrixFwd.h>
33+
#include <Core/GeometryPrimitives/Transform.h>
34+
#include <Core/Datatypes/Matrix.h>
35+
#include <Core/Datatypes/share.h>
36+
37+
namespace SCIRun {
38+
39+
namespace MatrixAlgorithms
40+
{
41+
#if SCIRUN4_TO_BE_ENABLED_LATER
42+
int cg_solve(const Matrix<double>& matrix, const ColumnMatrix& rhs, ColumnMatrix& lhs,
43+
double &err, int &niter,
44+
double max_error=1.e-6, int toomany=0,
45+
int useLhsAsGuess=0);
46+
int cg_solve(const Matrix<double>& matrix, const ColumnMatrix& rhs, ColumnMatrix& lhs);
47+
int cg_solve(const Matrix<double>& matrix, const DenseMatrix& rhs, DenseMatrix& lhs,
48+
double &err, int &niter,
49+
double max_error=1.e-6, int toomany=0,
50+
int useLhsAsGuess=0);
51+
int cg_solve(const Matrix<double>& matrix, const DenseMatrix& rhs, DenseMatrix& lhs);
52+
53+
int bicg_solve(const Matrix<double>& matrix, const ColumnMatrix& rhs, ColumnMatrix& lhs,
54+
double &err, int &niter,
55+
double max_error=1.e-6, int toomany=0,
56+
int useLhsAsGuess=0);
57+
int bicg_solve(const Matrix<double>& matrix, const ColumnMatrix& rhs, ColumnMatrix& lhs);
58+
int bicg_solve(const Matrix<double>& matrix, const DenseMatrix& rhs, DenseMatrix& lhs,
59+
double &err, int &niter,
60+
double max_error=1.e-6, int toomany=0,
61+
int useLhsAsGuess=0);
62+
int bicg_solve(const Matrix<double>& matrix, const DenseMatrix& rhs, DenseMatrix& lhs);
63+
#endif
64+
65+
SCISHARE Core::Geometry::Transform matrix_to_transform(const Core::Datatypes::Matrix& matrix);
66+
}
67+
}
68+
69+
#endif

0 commit comments

Comments
 (0)