1 /*************************************************************************
2 * *
3 * Vega FEM Simulation Library Version 3.1 *
4 * *
5 * "matrix" library , Copyright (C) 2007 CMU, 2009 MIT *
6 * All rights reserved. *
7 * *
8 * Code author: Jernej Barbic *
9 * http://www.jernejbarbic.com/code *
10 * Research: Jernej Barbic, Doug L. James, Jovan Popovic *
11 * Funding: NSF, Link Foundation, Singapore-MIT GAMBIT Game Lab *
12 * *
13 * This library is free software; you can redistribute it and/or *
14 * modify it under the terms of the BSD-style license that is *
15 * included with this library in the file LICENSE.txt *
16 * *
17 * This library is distributed in the hope that it will be useful, *
18 * but WITHOUT ANY WARRANTY; without even the implied warranty of *
19 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the file *
20 * LICENSE.TXT for more details. *
21 * *
22 *************************************************************************/
23
24 #include <stdio.h>
25 #include <stdlib.h>
26 #include <string.h>
27 #include <math.h>
28 #include "lapack-headers.h"
29 #include "matrixMacros.h"
30 #include "matrixBLAS.h"
31 #include "expokit_xgpadm.h"
32 #include "matrixExp.h"
33
34 #ifdef __APPLE__
35 #define DGPADM dgpadm_
36 #define SGPADM sgpadm_
37 #define INTEGER long int
38 #else
39 #define DGPADM dgpadm_
40 #define SGPADM sgpadm_
41 #define INTEGER int
42 #endif
43
44 template<bool C>
45 class _xgpadm {};
46
47 template<>
48 class _xgpadm<true> {
49 public:
f(INTEGER * ideg,INTEGER * m,float * t,float * H,INTEGER * ldh,float * wsp,INTEGER * lwsp,INTEGER * ipiv,INTEGER * iexph,INTEGER * ns,INTEGER * flag)50 inline static void f(INTEGER * ideg, INTEGER * m, float * t, float * H,
51 INTEGER * ldh, float * wsp, INTEGER * lwsp,
52 INTEGER * ipiv, INTEGER * iexph, INTEGER * ns, INTEGER * flag)
53 {
54 printf("Error: routine SGPADM is not supported by expokit.\n");
55 }
56 };
57
58
59 template<>
60 class _xgpadm<false> {
61 public:
f(INTEGER * ideg,INTEGER * m,double * t,double * H,INTEGER * ldh,double * wsp,INTEGER * lwsp,INTEGER * ipiv,INTEGER * iexph,INTEGER * ns,INTEGER * flag)62 inline static void f(INTEGER * ideg, INTEGER * m, double * t, double * H,
63 INTEGER * ldh, double * wsp, INTEGER * lwsp,
64 INTEGER * ipiv, INTEGER * iexph, INTEGER * ns, INTEGER * flag)
65 {
66 DGPADM(ideg, m, t, H, ldh, wsp, lwsp, ipiv, iexph, ns, flag);
67 }
68 };
69
70 // computes output = exp(mtx*t)
71 template<class real>
MatrixExp(int m,real * mtx,real t,real * output,int ideg_)72 int MatrixExp(int m, real * mtx, real t, real * output, int ideg_)
73 {
74 INTEGER ideg = ideg_;
75 INTEGER M = m;
76 INTEGER ldh = M;
77 INTEGER lwsp = 4*M*M+ideg+1;
78 real * wsp = (real*) malloc (sizeof(real) * lwsp);
79 INTEGER * ipiv = (INTEGER*) malloc (sizeof(INTEGER) * M);
80 INTEGER iexph = 0;
81 INTEGER ns = 0;
82 INTEGER flag = 0;
83
84 int mm = m*m;
85 real maxEntry = 0;
86 for(int i=0;i<mm; i++)
87 maxEntry = MAX(maxEntry, fabs(mtx[i]));
88
89 if (maxEntry != 0)
90 {
91 _xgpadm<sizeof(real)==sizeof(float)>::f(&ideg, &M, &t, mtx,
92 &ldh, wsp, &lwsp,
93 ipiv, &iexph, &ns, &flag);
94
95 //printf("iexph=%d ns=%d flag=%d\n", (int)iexph, (int)ns, (int)flag);
96
97 if (flag != 0)
98 {
99 printf("Error: xgpadm returned non-zero exit flag %d.\n", (int)flag);
100 return flag;
101 }
102
103 memcpy(output, &wsp[iexph-1], sizeof(real) * m * m);
104 }
105 else
106 {
107 // set result to identity matrix
108 memset(output, 0, sizeof(real) * m * m);
109 for(int i=0; i<m; i++)
110 output[ELT(m,i,i)] = 1.0;
111 }
112
113 free(wsp);
114
115 return 0;
116 }
117
118 // computes w = exp(mtx*t)*v
119 template<class real>
MatrixExpv(int m,real * mtx,real t,real * v,real * w,int ideg_)120 int MatrixExpv(int m, real * mtx, real t, real * v, real * w, int ideg_)
121 {
122 INTEGER ideg = ideg_;
123 INTEGER M = m;
124 INTEGER ldh = M;
125 INTEGER lwsp = 4*M*M+ideg+1;
126 real * wsp = (real*) malloc (sizeof(real) * lwsp);
127 INTEGER * ipiv = (INTEGER*) malloc (sizeof(INTEGER) * M);
128 INTEGER iexph = 0;
129 INTEGER ns = 0;
130 INTEGER flag = 0;
131
132 int mm = m*m;
133 real maxEntry = 0;
134 for(int i=0;i<mm; i++)
135 maxEntry = MAX(maxEntry, fabs(mtx[i]));
136
137 if (maxEntry != 0)
138 {
139 _xgpadm<sizeof(real)==sizeof(float)>::f(&ideg, &M, &t, mtx,
140 &ldh, wsp, &lwsp,
141 ipiv, &iexph, &ns, &flag);
142
143 //printf("iexph=%d ns=%d flag=%d\n", (int)iexph, (int)ns, (int)flag);
144
145 if (flag != 0)
146 {
147 printf("Error: xgpadm returned non-zero exit flag %d.\n", (int)flag);
148 return flag;
149 }
150
151 real * output = &wsp[iexph-1];
152 MultiplyMatrices(M, M, 1, output, v, w);
153 }
154 else
155 {
156 memcpy(w,v,sizeof(real)*M);
157 }
158
159 free(wsp);
160
161 return 0;
162 }
163
164 template int MatrixExpv(int m, float * mtx, float t, float * v, float * w, int ideg);
165 template int MatrixExpv(int m, double * mtx, double t, double * v, double * w, int ideg);
166
167 template int MatrixExp(int m, float * mtx, float t, float * output, int ideg_);
168 template int MatrixExp(int m, double * mtx, double t, double * output, int ideg_);
169