1 /* $Header: /var/cvs/mbdyn/mbdyn/mbdyn-1.0/libraries/libmbwrap/taucswrap.h,v 1.18 2017/01/12 14:44:25 masarati Exp $ */
2 /*
3  * HmFe (C) is a FEM analysis code.
4  *
5  * Copyright (C) 1996-2017
6  *
7  * Marco Morandini  <morandini@aero.polimi.it>
8  *
9  * Dipartimento di Ingegneria Aerospaziale - Politecnico di Milano
10  * via La Masa, 34 - 20156 Milano, Italy
11  * http://www.aero.polimi.it
12  *
13  * Changing this copyright notice is forbidden.
14  * This program is distributed in the hope that it will be useful,
15  * but WITHOUT ANY WARRANTY; without even the implied warranty of
16  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
17  *
18  */
19 /* December 2001
20  * Modified to add a Sparse matrix in row form and to implement methods
21  * to be used in the parallel MBDyn Solver.
22  *
23  * Copyright (C) 2001-2017
24  *
25  * Giuseppe Quaranta  <quaranta@aero.polimi.it>
26  *
27  * Dipartimento di Ingegneria Aerospaziale - Politecnico di Milano
28  * via La Masa, 34 - 20156 Milano, Italy
29  * http://www.aero.polimi.it
30  *
31  */
32 /*
33  * MBDyn (C) is a multibody analysis code.
34  * http://www.mbdyn.org
35  *
36  * Copyright (C) 1996-2017
37  *
38  * Pierangelo Masarati	<masarati@aero.polimi.it>
39  * Paolo Mantegazza	<mantegazza@aero.polimi.it>
40  *
41  * Dipartimento di Ingegneria Aerospaziale - Politecnico di Milano
42  * via La Masa, 34 - 20156 Milano, Italy
43  * http://www.aero.polimi.it
44  *
45  * Changing this copyright notice is forbidden.
46  *
47  * This program is free software; you can redistribute it and/or modify
48  * it under the terms of the GNU General Public License as published by
49  * the Free Software Foundation (version 2 of the License).
50  *
51  *
52  * This program is distributed in the hope that it will be useful,
53  * but WITHOUT ANY WARRANTY; without even the implied warranty of
54  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
55  * GNU General Public License for more details.
56  *
57  * You should have received a copy of the GNU General Public License
58  * along with this program; if not, write to the Free Software
59  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
60  */
61 
62 /*
63  * Taucs is used by permission; please read its Copyright,
64  * License and Availability note.
65  */
66 
67 #ifndef TaucsSparseSolutionManager_hh
68 #define TaucsSparseSolutionManager_hh
69 
70 #ifdef USE_TAUCS
71 
72 #include <iostream>
73 #include <vector>
74 
75 /*
76  * This is just to silence warnings.
77  * Do not use the taucs
78  * complex solver!!!!"
79 */
80 #include <complex>
81 
82 extern "C" {
83 #include <taucs.h>
84 }
85 
86 #include "myassert.h"
87 #include "mynewmem.h"
88 #include "ls.h"
89 #include "solman.h"
90 #include "spmapmh.h"
91 #include "ccmh.h"
92 
93 
94 /* TaucsSolver - begin */
95 
96 class TaucsSolver: public LinearSolver {
97 private:
98 	integer iSize;
99 	mutable taucs_ccs_matrix A;
100 	mutable doublereal *Axp;
101 	mutable integer *Aip;
102 	mutable integer *App;
103 	mutable bool Symbolic;
104 
105 	mutable void *Factorization;
106 
107 	void Factor(void);
108 
109 public:
110 	TaucsSolver(const integer &size);
111 	~TaucsSolver(void);
112 
113 	void Reset(void);
114 	void Solve(void) const;
115 
116 	void MakeCompactForm(SparseMatrixHandler&,
117 			std::vector<doublereal>& Ax,
118 			std::vector<integer>& Ar,
119 			std::vector<integer>& Ac,
120 			std::vector<integer>& Ap) const;
121 };
122 
123 /* TaucsSolver - end */
124 
125 /* TaucsSparseSolutionManager - begin */
126 
127 class TaucsSparseSolutionManager: public SolutionManager {
128 protected:
129 	mutable SpMapMatrixHandler A;
130 
131 	std::vector<doublereal> x;
132 	std::vector<doublereal> b;
133 
134 	mutable MyVectorHandler xVH, bVH;
135 
136 	std::vector<doublereal> Ax;
137 	std::vector<integer> Ai;
138 	std::vector<integer> Adummy;
139 	std::vector<integer> Ap;
140 
141 	/* Passa in forma di Compressed Column (callback per solve,
142 	 * richiesto da SpMap e CC Matrix Handler) */
143 	virtual void MakeCompressedColumnForm(void);
144 
145 	/* Backward Substitution */
146 	void BackSub(doublereal t_iniz = 0.);
147 
148 public:
149 	TaucsSparseSolutionManager(integer Dim);
150 	virtual ~TaucsSparseSolutionManager(void);
151 #ifdef DEBUG
IsValid(void)152 	virtual void IsValid(void) const {
153 		NO_OP;
154 	};
155 #endif /* DEBUG */
156 
157 	/* Inizializzatore generico */
158 	virtual void MatrReset(void);
159 
160 	/* Risolve il sistema Backward Substitution; fattorizza se necessario */
161 	virtual void Solve(void);
162 
163 	/* Rende disponibile l'handler per la matrice */
164 	virtual MatrixHandler* pMatHdl(void) const;
165 
166 	/* Rende disponibile l'handler per il termine noto */
167 	virtual MyVectorHandler* pResHdl(void) const;
168 
169 	/* Rende disponibile l'handler per la soluzione */
170 	virtual MyVectorHandler* pSolHdl(void) const;
171 };
172 
173 /* TaucsSparseSolutionManager - end */
174 
175 /* TaucsSparseCCSolutionManager - begin */
176 
177 template <class CC>
178 class TaucsSparseCCSolutionManager: public TaucsSparseSolutionManager {
179 protected:
180 	bool CCReady;
181 	CompactSparseMatrixHandler *Ac;
182 
183 	virtual void MatrReset(void);
184 	virtual void MakeCompressedColumnForm(void);
185 
186 public:
187 	TaucsSparseCCSolutionManager(integer Dim);
188 	virtual ~TaucsSparseCCSolutionManager(void);
189 
190 	/* Inizializzatore "speciale" */
191 	virtual void MatrInitialize(void);
192 
193 	/* Rende disponibile l'handler per la matrice */
194 	virtual MatrixHandler* pMatHdl(void) const;
195 };
196 
197 /* TaucsSparseCCSolutionManager - end */
198 
199 #endif /* USE_TAUCS */
200 
201 #endif /* TaucsSparseSolutionManager_hh */
202 
203