1 /*
2 XLiFE++ is an extended library of finite elements written in C++
3     Copyright (C) 2014  Lunéville, Eric; Kielbasiewicz, Nicolas; Lafranche, Yvon; Nguyen, Manh-Ha; Chambeyron, Colin
4 
5     This program is free software: you can redistribute it and/or modify
6     it under the terms of the GNU General Public License as published by
7     the Free Software Foundation, either version 3 of the License, or
8     (at your option) any later version.
9     This program is distributed in the hope that it will be useful,
10     but WITHOUT ANY WARRANTY; without even the implied warranty of
11     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
12     GNU General Public License for more details.
13     You should have received a copy of the GNU General Public License
14     along with this program.  If not, see <http://www.gnu.org/licenses/>.
15  */
16 
17 /*!
18   \file UmfPack.cpp
19   \authors Manh Ha NGUYEN, Y. Lafranche
20   \since 29 August 2013
21   \date 29 August 2013
22 
23   \brief Wrapper of some important functions of UMFPACK in order to solver (non)-symmetric problem
24 */
25 
26 #include "config.h"
27 #ifdef XLIFEPP_WITH_UMFPACK
28 
29 #include "UmfPack.hpp"
30 #include "UmfPackWrappers.hpp"
31 
32 namespace xlifepp
33 {
34 
35 //
36 // Specialization int/double
37 //
38 number_t UMFPACK<int, double>::infoSize = 90;
39 number_t UMFPACK<int, double>::controlSize = 30;
40 
41 
symbolic(int nRow,int nCol,const int Ap[],const int Ai[],const double Ax[],void ** Symbolic,const double Control[UMFPACK_CONTROL],double Info[UMFPACK_INFO])42 int UMFPACK<int, double>::symbolic(int nRow, int nCol,
43                     const int Ap[], const int Ai[],
44                     const double Ax[], void** Symbolic,
45                     const double Control[UMFPACK_CONTROL], double Info[UMFPACK_INFO])
46 {
47     return (DISYMBOLIC(nRow, nCol, Ap, Ai, Ax, Symbolic, Control, Info));
48 }
49 
numeric(const int Ap[],const int Ai[],const double Ax[],void * Symbolic,void ** Numeric,const double Control[UMFPACK_CONTROL],double Info[UMFPACK_INFO])50 int UMFPACK<int, double>::numeric(const int Ap[], const int Ai[],
51                     const double Ax[], void* Symbolic, void** Numeric,
52                     const double Control[UMFPACK_CONTROL], double Info[UMFPACK_INFO])
53 {
54     return (DINUMERIC(Ap, Ai, Ax, Symbolic, Numeric, Control, Info));
55 }
56 
solve(int sys,const int Ap[],const int Ai[],const double Ax[],double X[],const double B[],void * Numeric,const double Control[UMFPACK_CONTROL],double Info[UMFPACK_INFO])57 int UMFPACK<int, double>::solve(int sys, const int Ap[], const int Ai[],
58         const double Ax[], double X[], const double B[], void* Numeric,
59         const double Control[UMFPACK_CONTROL], double Info[UMFPACK_INFO])
60 {
61     return (DISOLVE(sys, Ap, Ai, Ax, X, B, Numeric, Control, Info));
62 }
63 
freeSymbolic(void ** Symbolic)64 void UMFPACK<int, double>::freeSymbolic(void** Symbolic)
65 {
66     DIFREESYMBOLIC(Symbolic);
67 }
68 
freeNumeric(void ** Numeric)69 void UMFPACK<int, double>::freeNumeric(void** Numeric)
70 {
71     DIFREENUMERIC(Numeric);
72 }
73 
getLunz(int * lnz,int * unz,int * nRow,int * nCol,int * nzUdiag,void * Numeric)74 int UMFPACK<int, double>::getLunz(int* lnz, int* unz, int* nRow, int* nCol, int* nzUdiag, void* Numeric)
75 {
76     return DIGETLUNZ(lnz, unz, nRow, nCol, nzUdiag, Numeric);
77 }
78 
getNumeric(int Lp[],int Lj[],double Lx[],int Up[],int Ui[],double Ux[],int P[],int Q[],double Dx[],int * doRecip,double Rs[],void * Numeric)79 int UMFPACK<int, double>::getNumeric(int Lp[], int Lj[], double Lx[], int Up[], int Ui[], double Ux[],
80         int P[], int Q[], double Dx[], int* doRecip, double Rs[], void* Numeric)
81 {
82     return DIGETNUMERIC(Lp, Lj, Lx, Up, Ui, Ux, P, Q, Dx, doRecip, Rs, Numeric);
83 }
84 
getDeterminant(double * Mx,double * Ex,void * NumericHandle,double UserInfo[UMFPACK_INFO])85 int UMFPACK<int, double>::getDeterminant(double* Mx, double* Ex, void* NumericHandle, double UserInfo[UMFPACK_INFO])
86 {
87     return DIGETDET(Mx, Ex, NumericHandle, UserInfo);
88 }
89 
90 //
91 // Specialization int/complex<double>
92 //
93 
symbolic(int nRow,int nCol,const int Ap[],const int Ai[],const std::complex<double> Az[],void ** Symbolic,const double Control[UMFPACK_CONTROL],double Info[UMFPACK_INFO])94 int UMFPACK<int, std::complex<double> >::symbolic(int nRow, int nCol,
95         const int Ap[], const int Ai[],
96         const std::complex<double> Az[], void** Symbolic,
97         const double Control[UMFPACK_CONTROL], double Info[UMFPACK_INFO])
98 {
99     double* null = (double*)NULL;
100     const double* Axy=reinterpret_cast<const double*>(Az); //packed complex
101     return (ZISYMBOLIC(nRow, nCol, Ap, Ai, Axy, null, Symbolic, Control, Info));
102 }
103 
numeric(const int Ap[],const int Ai[],const std::complex<double> Az[],void * Symbolic,void ** Numeric,const double Control[UMFPACK_CONTROL],double Info[UMFPACK_INFO])104 int UMFPACK<int, std::complex<double> >::numeric(const int Ap[], const int Ai[],
105         const std::complex<double> Az[], void* Symbolic, void** Numeric,
106         const double Control[UMFPACK_CONTROL], double Info[UMFPACK_INFO])
107 {
108     double* null = (double*)NULL;
109     const double* Axy=reinterpret_cast<const double*>(Az);  //packed complex
110     return ZINUMERIC(Ap, Ai, Axy, null, Symbolic, Numeric, Control, Info);
111 }
112 
solve(int sys,const int Ap[],const int Ai[],const std::complex<double> Az[],std::complex<double> X[],const std::complex<double> B[],void * Numeric,const double Control[UMFPACK_CONTROL],double Info[UMFPACK_INFO])113 int UMFPACK<int, std::complex<double> >::solve(int sys, const int Ap[], const int Ai[],
114         const std::complex<double> Az[], std::complex<double> X[], const std::complex<double> B[],
115         void* Numeric, const double Control[UMFPACK_CONTROL], double Info[UMFPACK_INFO])
116 {
117     double* null = (double*)NULL;
118     const double* Axy=reinterpret_cast<const double*>(Az);  //packed complex
119     const double* Bxy=reinterpret_cast<const double*>(B);   //packed complex
120     double* Xxy=reinterpret_cast<double*>(X);               //packed complex
121     return (ZISOLVE(sys, Ap, Ai, Axy, null, Xxy, null, Bxy, null, Numeric, Control, Info));
122 }
123 
solveCR(int sys,const int Ap[],const int Ai[],const std::complex<double> Az[],std::complex<double> X[],const double Bx[],const double Bz[],void * Numeric,const double Control[UMFPACK_CONTROL],double Info[UMFPACK_INFO])124 int UMFPACK<int, std::complex<double> >::solveCR(int sys, const int Ap[], const int Ai[],
125         const std::complex<double> Az[], std::complex<double> X[], const double Bx[], const double Bz[],
126         void* Numeric, const double Control[UMFPACK_CONTROL], double Info[UMFPACK_INFO])
127 {
128     double* null = (double*)NULL;
129     const double* Axy=reinterpret_cast<const double*>(Az);  //packed complex
130     double* Xxy=reinterpret_cast<double*>(X);               //packed complex
131     return (ZISOLVE(sys, Ap, Ai, Axy, null, Xxy, null, Bx, Bz, Numeric, Control, Info));
132 }
133 
freeSymbolic(void ** Symbolic)134 void UMFPACK<int, std::complex<double> >::freeSymbolic(void** Symbolic)
135 {
136     ZIFREESYMBOLIC(Symbolic);
137 }
138 
freeNumeric(void ** Numeric)139 void UMFPACK<int, std::complex<double> >::freeNumeric(void** Numeric)
140 {
141     ZIFREENUMERIC(Numeric);
142 }
143 
getLunz(int * lnz,int * unz,int * nRow,int * nCol,int * nzUdiag,void * Numeric)144 int UMFPACK<int, std::complex<double> >::getLunz(int* lnz, int* unz, int* nRow, int* nCol, int* nzUdiag, void* Numeric)
145 {
146     return ZIGETLUNZ(lnz, unz, nRow, nCol, nzUdiag, Numeric);
147 }
148 
getNumeric(int Lp[],int Lj[],std::complex<double> Lz[],int Up[],int Ui[],std::complex<double> Uz[],int P[],int Q[],std::complex<double> Dz[],int * doRecip,double Rs[],void * Numeric)149 int UMFPACK<int, std::complex<double> >::getNumeric(int Lp[], int Lj[], std::complex<double> Lz[], int Up[], int Ui[], std::complex<double> Uz[],
150         int P[], int Q[], std::complex<double> Dz[], int *doRecip, double Rs[], void* Numeric)
151 {
152     double* null = (double*)NULL;
153     double* Lxy=reinterpret_cast<double*>(Lz);  //packed complex
154     double* Uxy=reinterpret_cast<double*>(Uz);  //packed complex
155     double* Dxy=reinterpret_cast<double*>(Dz);  //packed complex
156     return ZIGETNUMERIC(Lp, Lj, Lxy, null, Up, Ui, Uxy, null, P, Q, Dxy, null, doRecip, Rs, Numeric);
157 }
158 
getDeterminant(std::complex<double> * Mx,double * Ex,void * NumericHandle,double UserInfo[UMFPACK_INFO])159 int UMFPACK<int, std::complex<double> >::getDeterminant(std::complex<double>* Mx, double* Ex, void* NumericHandle, double UserInfo[UMFPACK_INFO])
160 {
161     double* null = (double*)NULL;
162     return ZIGETDET(reinterpret_cast<double*>(Mx), null, Ex, NumericHandle, UserInfo);
163 }
164 
165 //
166 // Specialization long long/double
167 //
168 
symbolic(SuiteSparseLong nRow,SuiteSparseLong nCol,const SuiteSparseLong Ap[],const SuiteSparseLong Ai[],const double Ax[],void ** Symbolic,const double Control[UMFPACK_CONTROL],double Info[UMFPACK_INFO])169 SuiteSparseLong UMFPACK<SuiteSparseLong, double>::symbolic(SuiteSparseLong nRow, SuiteSparseLong nCol,
170                     const SuiteSparseLong Ap[], const SuiteSparseLong Ai[],
171                     const double Ax[], void** Symbolic,
172                     const double Control[UMFPACK_CONTROL], double Info[UMFPACK_INFO])
173 {
174     return (DLSYMBOLIC(nRow, nCol, Ap, Ai, Ax, Symbolic, Control, Info));
175 }
176 
numeric(const SuiteSparseLong Ap[],const SuiteSparseLong Ai[],const double Ax[],void * Symbolic,void ** Numeric,const double Control[UMFPACK_CONTROL],double Info[UMFPACK_INFO])177 SuiteSparseLong UMFPACK<SuiteSparseLong, double>::numeric(const SuiteSparseLong Ap[], const SuiteSparseLong Ai[],
178                     const double Ax[], void* Symbolic, void** Numeric,
179                     const double Control[UMFPACK_CONTROL], double Info[UMFPACK_INFO])
180 {
181     return (DLNUMERIC(Ap, Ai, Ax, Symbolic, Numeric, Control, Info));
182 }
183 
solve(SuiteSparseLong sys,const SuiteSparseLong Ap[],const SuiteSparseLong Ai[],const double Ax[],double X[],const double B[],void * Numeric,const double Control[UMFPACK_CONTROL],double Info[UMFPACK_INFO])184 SuiteSparseLong UMFPACK<SuiteSparseLong, double>::solve(SuiteSparseLong sys, const SuiteSparseLong Ap[], const SuiteSparseLong Ai[],
185         const double Ax[], double X[], const double B[], void* Numeric,
186         const double Control[UMFPACK_CONTROL], double Info[UMFPACK_INFO])
187 {
188     return (DLSOLVE(sys, Ap, Ai, Ax, X, B, Numeric, Control, Info));
189 }
190 
freeSymbolic(void ** Symbolic)191 void UMFPACK<SuiteSparseLong, double>::freeSymbolic(void** Symbolic)
192 {
193     DLFREESYMBOLIC(Symbolic);
194 }
195 
freeNumeric(void ** Numeric)196 void UMFPACK<SuiteSparseLong, double>::freeNumeric(void** Numeric)
197 {
198     DLFREENUMERIC(Numeric);
199 }
200 
getLunz(SuiteSparseLong * lnz,SuiteSparseLong * unz,SuiteSparseLong * nRow,SuiteSparseLong * nCol,SuiteSparseLong * nzUdiag,void * Numeric)201 SuiteSparseLong UMFPACK<SuiteSparseLong, double>::getLunz(SuiteSparseLong* lnz, SuiteSparseLong* unz, SuiteSparseLong* nRow, SuiteSparseLong* nCol, SuiteSparseLong* nzUdiag, void* Numeric)
202 {
203     return DLGETLUNZ(lnz, unz, nRow, nCol, nzUdiag, Numeric);
204 }
205 
getNumeric(SuiteSparseLong Lp[],SuiteSparseLong Lj[],double Lx[],SuiteSparseLong Up[],SuiteSparseLong Ui[],double Ux[],SuiteSparseLong P[],SuiteSparseLong Q[],double Dx[],SuiteSparseLong * doRecip,double Rs[],void * Numeric)206 SuiteSparseLong UMFPACK<SuiteSparseLong, double>::getNumeric(SuiteSparseLong Lp[], SuiteSparseLong Lj[], double Lx[], SuiteSparseLong Up[], SuiteSparseLong Ui[], double Ux[],
207         SuiteSparseLong P[], SuiteSparseLong Q[], double Dx[], SuiteSparseLong* doRecip, double Rs[], void* Numeric)
208 {
209     return DLGETNUMERIC(Lp, Lj, Lx, Up, Ui, Ux, P, Q, Dx, doRecip, Rs, Numeric);
210 }
211 
getDeterminant(double * Mx,double * Ex,void * NumericHandle,double UserInfo[UMFPACK_INFO])212 SuiteSparseLong UMFPACK<SuiteSparseLong, double>::getDeterminant(double* Mx, double* Ex, void* NumericHandle, double UserInfo[UMFPACK_INFO])
213 {
214     return DLGETDET(Mx, Ex, NumericHandle, UserInfo);
215 }
216 
217 //
218 // Specialization long long/complex<double>
219 //
220 
symbolic(SuiteSparseLong nRow,SuiteSparseLong nCol,const SuiteSparseLong Ap[],const SuiteSparseLong Ai[],const std::complex<double> Az[],void ** Symbolic,const double Control[UMFPACK_CONTROL],double Info[UMFPACK_INFO])221 SuiteSparseLong UMFPACK<SuiteSparseLong, std::complex<double> >::symbolic(SuiteSparseLong nRow, SuiteSparseLong nCol,
222         const SuiteSparseLong Ap[], const SuiteSparseLong Ai[],
223         const std::complex<double> Az[], void** Symbolic,
224         const double Control[UMFPACK_CONTROL], double Info[UMFPACK_INFO])
225 {
226     double* null = (double*)NULL;
227     const double* Axy=reinterpret_cast<const double*>(Az); //packed complex
228     return (ZLSYMBOLIC(nRow, nCol, Ap, Ai, Axy, null, Symbolic, Control, Info));
229 }
230 
numeric(const SuiteSparseLong Ap[],const SuiteSparseLong Ai[],const std::complex<double> Az[],void * Symbolic,void ** Numeric,const double Control[UMFPACK_CONTROL],double Info[UMFPACK_INFO])231 SuiteSparseLong UMFPACK<SuiteSparseLong, std::complex<double> >::numeric(const SuiteSparseLong Ap[], const SuiteSparseLong Ai[],
232         const std::complex<double> Az[], void* Symbolic, void** Numeric,
233         const double Control[UMFPACK_CONTROL], double Info[UMFPACK_INFO])
234 {
235     double* null = (double*)NULL;
236     const double* Axy=reinterpret_cast<const double*>(Az);  //packed complex
237     return ZLNUMERIC(Ap, Ai, Axy, null, Symbolic, Numeric, Control, Info);
238 }
239 
solve(SuiteSparseLong sys,const SuiteSparseLong Ap[],const SuiteSparseLong Ai[],const std::complex<double> Az[],std::complex<double> X[],const std::complex<double> B[],void * Numeric,const double Control[UMFPACK_CONTROL],double Info[UMFPACK_INFO])240 SuiteSparseLong UMFPACK<SuiteSparseLong, std::complex<double> >::solve(SuiteSparseLong sys, const SuiteSparseLong Ap[], const SuiteSparseLong Ai[],
241         const std::complex<double> Az[], std::complex<double> X[], const std::complex<double> B[],
242         void* Numeric, const double Control[UMFPACK_CONTROL], double Info[UMFPACK_INFO])
243 {
244     double* null = (double*)NULL;
245     const double* Axy=reinterpret_cast<const double*>(Az);  //packed complex
246     const double* Bxy=reinterpret_cast<const double*>(B);   //packed complex
247     double* Xxy=reinterpret_cast<double*>(X);               //packed complex
248     return (ZLSOLVE(sys, Ap, Ai, Axy, null, Xxy, null, Bxy, null, Numeric, Control, Info));
249 }
250 
solveCR(SuiteSparseLong sys,const SuiteSparseLong Ap[],const SuiteSparseLong Ai[],const std::complex<double> Az[],std::complex<double> X[],const double Bx[],const double Bz[],void * Numeric,const double Control[UMFPACK_CONTROL],double Info[UMFPACK_INFO])251 SuiteSparseLong UMFPACK<SuiteSparseLong, std::complex<double> >::solveCR(SuiteSparseLong sys, const SuiteSparseLong Ap[], const SuiteSparseLong Ai[],
252         const std::complex<double> Az[], std::complex<double> X[], const double Bx[], const double Bz[],
253         void* Numeric, const double Control[UMFPACK_CONTROL], double Info[UMFPACK_INFO])
254 {
255     double* null = (double*)NULL;
256     const double* Axy=reinterpret_cast<const double*>(Az);  //packed complex
257     double* Xxy=reinterpret_cast<double*>(X);               //packed complex
258     return (ZLSOLVE(sys, Ap, Ai, Axy, null, Xxy, null, Bx, Bz, Numeric, Control, Info));
259 }
260 
freeSymbolic(void ** Symbolic)261 void UMFPACK<SuiteSparseLong, std::complex<double> >::freeSymbolic(void** Symbolic)
262 {
263     ZLFREESYMBOLIC(Symbolic);
264 }
265 
freeNumeric(void ** Numeric)266 void UMFPACK<SuiteSparseLong, std::complex<double> >::freeNumeric(void** Numeric)
267 {
268     ZLFREENUMERIC(Numeric);
269 }
270 
getLunz(SuiteSparseLong * lnz,SuiteSparseLong * unz,SuiteSparseLong * nRow,SuiteSparseLong * nCol,SuiteSparseLong * nzUdiag,void * Numeric)271 SuiteSparseLong UMFPACK<SuiteSparseLong, std::complex<double> >::getLunz(SuiteSparseLong* lnz, SuiteSparseLong* unz, SuiteSparseLong* nRow, SuiteSparseLong* nCol, SuiteSparseLong* nzUdiag, void* Numeric)
272 {
273     return ZLGETLUNZ(lnz, unz, nRow, nCol, nzUdiag, Numeric);
274 }
275 
getNumeric(SuiteSparseLong Lp[],SuiteSparseLong Lj[],std::complex<double> Lz[],SuiteSparseLong Up[],SuiteSparseLong Ui[],std::complex<double> Uz[],SuiteSparseLong P[],SuiteSparseLong Q[],std::complex<double> Dz[],SuiteSparseLong * doRecip,double Rs[],void * Numeric)276 SuiteSparseLong UMFPACK<SuiteSparseLong, std::complex<double> >::getNumeric(SuiteSparseLong Lp[], SuiteSparseLong Lj[], std::complex<double> Lz[], SuiteSparseLong Up[], SuiteSparseLong Ui[], std::complex<double> Uz[],
277         SuiteSparseLong P[], SuiteSparseLong Q[], std::complex<double> Dz[], SuiteSparseLong *doRecip, double Rs[], void* Numeric)
278 {
279     double* null = (double*)NULL;
280     double* Lxy=reinterpret_cast<double*>(Lz);  //packed complex
281     double* Uxy=reinterpret_cast<double*>(Uz);  //packed complex
282     double* Dxy=reinterpret_cast<double*>(Dz);  //packed complex
283     return ZLGETNUMERIC(Lp, Lj, Lxy, null, Up, Ui, Uxy, null, P, Q, Dxy, null, doRecip, Rs, Numeric);
284 }
285 
getDeterminant(std::complex<double> * Mx,double * Ex,void * NumericHandle,double UserInfo[UMFPACK_INFO])286 SuiteSparseLong UMFPACK<SuiteSparseLong, std::complex<double> >::getDeterminant(std::complex<double>* Mx, double* Ex, void* NumericHandle, double UserInfo[UMFPACK_INFO])
287 {
288     double* null = (double*)NULL;
289     return ZLGETDET(reinterpret_cast<double*>(Mx), null, Ex, NumericHandle, UserInfo);
290 }
291 
292 
293 }
294 #endif  //XLIFEPP_WITH_UMFPACK
295 
296