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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2742 - (hide annotations)
Thu Nov 12 06:03:37 2009 UTC (10 years ago) by jfenwick
File size: 2134 byte(s)
Merging changes from the lapack branch.

The inverse() operation has been moved into c++. [No lazy support for this operation yet.]
Optional Lapack support has been added for matrices larger than 3x3. 
service0 is set to use mkl_lapack.



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

  ViewVC Help
Powered by ViewVC 1.1.26