/[escript]/trunk/escript/src/LapackInverseHelper.cpp
ViewVC logotype

Contents of /trunk/escript/src/LapackInverseHelper.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4369 - (show annotations)
Fri Apr 19 02:32:34 2013 UTC (5 years, 9 months ago) by jfenwick
File size: 2259 byte(s)
fix problems revealed on freebsd
1
2 /*****************************************************************************
3 *
4 * Copyright (c) 2009-2013 by University of Queensland
5 * http://www.uq.edu.au
6 *
7 * Primary Business: Queensland, Australia
8 * Licensed under the Open Software License version 3.0
9 * http://www.opensource.org/licenses/osl-3.0.php
10 *
11 * Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12 * Development since 2012 by School of Earth Sciences
13 *
14 *****************************************************************************/
15
16 #include "LapackInverseHelper.h"
17
18 #ifdef USE_LAPACK
19
20 #ifdef MKL_LAPACK
21 #include <mkl_lapack.h>
22 #else // assuming clapack
23 extern "C"
24 {
25 #include <clapack.h>
26 }
27 #endif
28
29 #endif
30
31
32 using namespace escript;
33
34 #define SUCCESS 0
35 #define NEEDLAPACK 5
36 #define ERRFACTORISE 6
37 #define ERRINVERT 7
38
39 LapackInverseHelper::LapackInverseHelper(int N)
40 {
41 piv=0;
42 work=0;
43 lwork=0;
44 this->N=N;
45 #ifdef USE_LAPACK
46 piv=new int[N];
47 int blocksize=64; // this is arbitrary. For implementations that require work array
48 // maybe we should look into the Lapack ILAENV function
49 #ifdef MKL_LAPACK
50 int minus1=-1;
51 double dummyd=0;
52 int result=0;
53 dgetri(&N, &dummyd, &N, &N, &dummyd, &minus1, &result); // The only param that matters is the second -1
54 if (result==0)
55 {
56 blocksize=static_cast<int>(dummyd);
57 }
58 // If there is an error, then fail silently.
59 // Why: I need this to be threadsafe so I can't throw and If this call fails,
60 // It will fail again when we try to invert the matrix
61 #endif
62 lwork=N*blocksize;
63 work=new double[lwork];
64 #endif
65 }
66
67 LapackInverseHelper::~LapackInverseHelper()
68 {
69 if (piv!=0)
70 {
71 delete[] piv;
72 }
73 if (work!=0)
74 {
75 delete[] work;
76 }
77 }
78
79 int
80 LapackInverseHelper::invert(double* matrix)
81 {
82 #ifndef USE_LAPACK
83 return NEEDLAPACK;
84 #else
85 #ifdef MKL_LAPACK
86 int res=0;
87 int size=N;
88 dgetrf(&N,&N,matrix,&N,piv,&res);
89 if (res!=0)
90 {
91 return ERRFACTORISE;
92 }
93 dgetri(&N, matrix, &N, piv, work, &lwork, &res);
94 if (res!=0)
95 {
96 return ERRINVERT;
97 }
98 #else // assuming clapack
99 int res=clapack_dgetrf(CblasColMajor, N,N,matrix,N, piv);
100 if (res!=0)
101 {
102 return ERRFACTORISE;
103 }
104 res=clapack_dgetri(CblasColMajor ,N,matrix,N , piv);
105 if (res!=0)
106 {
107 return ERRINVERT;
108 }
109 #endif
110 return SUCCESS;
111 #endif
112 }
113

  ViewVC Help
Powered by ViewVC 1.1.26