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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4369 - (hide annotations)
Fri Apr 19 02:32:34 2013 UTC (6 years ago) by jfenwick
File size: 2259 byte(s)
fix problems revealed on freebsd
1 jfenwick 2742
2 jfenwick 3981 /*****************************************************************************
3 jfenwick 2742 *
4 jfenwick 4154 * Copyright (c) 2009-2013 by University of Queensland
5 jfenwick 3981 * http://www.uq.edu.au
6 jfenwick 2742 *
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 jfenwick 3981 * Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12     * Development since 2012 by School of Earth Sciences
13     *
14     *****************************************************************************/
15 jfenwick 2742
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 jfenwick 4369

  ViewVC Help
Powered by ViewVC 1.1.26