/[escript]/branches/subworld2/escriptcore/src/LapackInverseHelper.cpp
ViewVC logotype

Contents of /branches/subworld2/escriptcore/src/LapackInverseHelper.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2881 - (show annotations)
Thu Jan 28 02:03:15 2010 UTC (9 years, 3 months ago) by jfenwick
Original Path: trunk/escript/src/LapackInverseHelper.cpp
File size: 2139 byte(s)
Don't panic.
Updating copyright stamps

1
2 /*******************************************************
3 *
4 * Copyright (c) 2009-2010 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