Actual source code: fnexp.c
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-, Universitat Politecnica de Valencia, Spain
6: This file is part of SLEPc.
7: SLEPc is distributed under a 2-clause BSD license (see LICENSE).
8: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
9: */
10: /*
11: Exponential function exp(x)
12: */
14: #include <slepc/private/fnimpl.h>
15: #include <slepcblaslapack.h>
17: static PetscErrorCode FNEvaluateFunction_Exp(FN fn,PetscScalar x,PetscScalar *y)
18: {
19: PetscFunctionBegin;
20: *y = PetscExpScalar(x);
21: PetscFunctionReturn(PETSC_SUCCESS);
22: }
24: static PetscErrorCode FNEvaluateDerivative_Exp(FN fn,PetscScalar x,PetscScalar *y)
25: {
26: PetscFunctionBegin;
27: *y = PetscExpScalar(x);
28: PetscFunctionReturn(PETSC_SUCCESS);
29: }
31: #define MAX_PADE 6
33: static PetscErrorCode FNEvaluateFunctionMat_Exp_Pade(FN fn,Mat A,Mat B)
34: {
35: PetscBLASInt n=0,ld,ld2,*ipiv,inc=1;
36: PetscInt m,j,k,sexp;
37: PetscBool odd;
38: const PetscInt p=MAX_PADE;
39: PetscReal c[MAX_PADE+1],s,*rwork;
40: PetscScalar scale,mone=-1.0,one=1.0,two=2.0,zero=0.0;
41: PetscScalar *Ba,*As,*A2,*Q,*P,*W,*aux;
42: const PetscScalar *Aa;
44: PetscFunctionBegin;
45: PetscCall(MatDenseGetArrayRead(A,&Aa));
46: PetscCall(MatDenseGetArray(B,&Ba));
47: PetscCall(MatGetSize(A,&m,NULL));
48: PetscCall(PetscBLASIntCast(m,&n));
49: ld = n;
50: ld2 = ld*ld;
51: P = Ba;
52: PetscCall(PetscMalloc6(m*m,&Q,m*m,&W,m*m,&As,m*m,&A2,ld,&rwork,ld,&ipiv));
53: PetscCall(PetscArraycpy(As,Aa,ld2));
55: /* Pade' coefficients */
56: c[0] = 1.0;
57: for (k=1;k<=p;k++) c[k] = c[k-1]*(p+1-k)/(k*(2*p+1-k));
59: /* Scaling */
60: s = LAPACKlange_("I",&n,&n,As,&ld,rwork);
61: PetscCall(PetscLogFlops(1.0*n*n));
62: if (s>0.5) {
63: sexp = PetscMax(0,(int)(PetscLogReal(s)/PetscLogReal(2.0))+2);
64: scale = PetscPowRealInt(2.0,-sexp);
65: PetscCallBLAS("BLASscal",BLASscal_(&ld2,&scale,As,&inc));
66: PetscCall(PetscLogFlops(1.0*n*n));
67: } else sexp = 0;
69: /* Horner evaluation */
70: PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&n,&n,&n,&one,As,&ld,As,&ld,&zero,A2,&ld));
71: PetscCall(PetscLogFlops(2.0*n*n*n));
72: PetscCall(PetscArrayzero(Q,ld2));
73: PetscCall(PetscArrayzero(P,ld2));
74: for (j=0;j<n;j++) {
75: Q[j+j*ld] = c[p];
76: P[j+j*ld] = c[p-1];
77: }
79: odd = PETSC_TRUE;
80: for (k=p-1;k>0;k--) {
81: if (odd) {
82: PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&n,&n,&n,&one,Q,&ld,A2,&ld,&zero,W,&ld));
83: SlepcSwap(Q,W,aux);
84: for (j=0;j<n;j++) Q[j+j*ld] += c[k-1];
85: odd = PETSC_FALSE;
86: } else {
87: PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&n,&n,&n,&one,P,&ld,A2,&ld,&zero,W,&ld));
88: SlepcSwap(P,W,aux);
89: for (j=0;j<n;j++) P[j+j*ld] += c[k-1];
90: odd = PETSC_TRUE;
91: }
92: PetscCall(PetscLogFlops(2.0*n*n*n));
93: }
94: /*if (odd) {
95: PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&n,&n,&n,&one,Q,&ld,As,&ld,&zero,W,&ld));
96: SlepcSwap(Q,W,aux);
97: PetscCallBLAS("BLASaxpy",BLASaxpy_(&ld2,&mone,P,&inc,Q,&inc));
98: PetscCallLAPACKInfo("LAPACKgesv",LAPACKgesv_(&n,&n,Q,&ld,ipiv,P,&ld,&info));
99: PetscCallBLAS("BLASscal",BLASscal_(&ld2,&two,P,&inc));
100: for (j=0;j<n;j++) P[j+j*ld] += 1.0;
101: PetscCallBLAS("BLASscal",BLASscal_(&ld2,&mone,P,&inc));
102: } else {*/
103: PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&n,&n,&n,&one,P,&ld,As,&ld,&zero,W,&ld));
104: SlepcSwap(P,W,aux);
105: PetscCallBLAS("BLASaxpy",BLASaxpy_(&ld2,&mone,P,&inc,Q,&inc));
106: PetscCallLAPACKInfo("LAPACKgesv",LAPACKgesv_(&n,&n,Q,&ld,ipiv,P,&ld,&info));
107: PetscCallBLAS("BLASscal",BLASscal_(&ld2,&two,P,&inc));
108: for (j=0;j<n;j++) P[j+j*ld] += 1.0;
109: /*}*/
110: PetscCall(PetscLogFlops(2.0*n*n*n+2.0*n*n*n/3.0+4.0*n*n));
112: for (k=1;k<=sexp;k++) {
113: PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&n,&n,&n,&one,P,&ld,P,&ld,&zero,W,&ld));
114: PetscCall(PetscArraycpy(P,W,ld2));
115: }
116: if (P!=Ba) PetscCall(PetscArraycpy(Ba,P,ld2));
117: PetscCall(PetscLogFlops(2.0*n*n*n*sexp));
119: PetscCall(PetscFree6(Q,W,As,A2,rwork,ipiv));
120: PetscCall(MatDenseRestoreArrayRead(A,&Aa));
121: PetscCall(MatDenseRestoreArray(B,&Ba));
122: PetscFunctionReturn(PETSC_SUCCESS);
123: }
125: #if PetscDefined(HAVE_COMPLEX)
126: /*
127: * Set scaling factor (s) and Pade degree (k,m)
128: */
129: static PetscErrorCode sexpm_params(PetscReal nrm,PetscInt *s,PetscInt *k,PetscInt *m)
130: {
131: PetscFunctionBegin;
132: if (nrm>1) {
133: if (nrm<200) {*s = 4; *k = 5; *m = *k-1;}
134: else if (nrm<1e4) {*s = 4; *k = 4; *m = *k+1;}
135: else if (nrm<1e6) {*s = 4; *k = 3; *m = *k+1;}
136: else if (nrm<1e9) {*s = 3; *k = 3; *m = *k+1;}
137: else if (nrm<1e11) {*s = 2; *k = 3; *m = *k+1;}
138: else if (nrm<1e12) {*s = 2; *k = 2; *m = *k+1;}
139: else if (nrm<1e14) {*s = 2; *k = 1; *m = *k+1;}
140: else {*s = 1; *k = 1; *m = *k+1;}
141: } else { /* nrm<1 */
142: if (nrm>0.5) {*s = 4; *k = 4; *m = *k-1;}
143: else if (nrm>0.3) {*s = 3; *k = 4; *m = *k-1;}
144: else if (nrm>0.15) {*s = 2; *k = 4; *m = *k-1;}
145: else if (nrm>0.07) {*s = 1; *k = 4; *m = *k-1;}
146: else if (nrm>0.01) {*s = 0; *k = 4; *m = *k-1;}
147: else if (nrm>3e-4) {*s = 0; *k = 3; *m = *k-1;}
148: else if (nrm>1e-5) {*s = 0; *k = 3; *m = 0;}
149: else if (nrm>1e-8) {*s = 0; *k = 2; *m = 0;}
150: else {*s = 0; *k = 1; *m = 0;}
151: }
152: PetscFunctionReturn(PETSC_SUCCESS);
153: }
155: /*
156: * Partial fraction form coefficients.
157: * If query, the function returns the size necessary to store the coefficients.
158: */
159: static PetscErrorCode getcoeffs(PetscInt k,PetscInt m,PetscComplex *r,PetscComplex *q,PetscComplex *remain,PetscBool query)
160: {
161: PetscInt i;
162: const PetscComplex /* m == k+1 */
163: p1r4[5] = {-1.582680186458572e+01 - 2.412564578224361e+01*PETSC_i,
164: -1.582680186458572e+01 + 2.412564578224361e+01*PETSC_i,
165: 1.499984465975511e+02 + 6.804227952202417e+01*PETSC_i,
166: 1.499984465975511e+02 - 6.804227952202417e+01*PETSC_i,
167: -2.733432894659307e+02 },
168: p1q4[5] = { 3.655694325463550e+00 + 6.543736899360086e+00*PETSC_i,
169: 3.655694325463550e+00 - 6.543736899360086e+00*PETSC_i,
170: 5.700953298671832e+00 + 3.210265600308496e+00*PETSC_i,
171: 5.700953298671832e+00 - 3.210265600308496e+00*PETSC_i,
172: 6.286704751729261e+00 },
173: p1r3[4] = {-1.130153999597152e+01 + 1.247167585025031e+01*PETSC_i,
174: -1.130153999597152e+01 - 1.247167585025031e+01*PETSC_i,
175: 1.330153999597152e+01 - 6.007173273704750e+01*PETSC_i,
176: 1.330153999597152e+01 + 6.007173273704750e+01*PETSC_i},
177: p1q3[4] = { 3.212806896871536e+00 + 4.773087433276636e+00*PETSC_i,
178: 3.212806896871536e+00 - 4.773087433276636e+00*PETSC_i,
179: 4.787193103128464e+00 + 1.567476416895212e+00*PETSC_i,
180: 4.787193103128464e+00 - 1.567476416895212e+00*PETSC_i},
181: p1r2[3] = { 7.648749087422928e+00 + 4.171640244747463e+00*PETSC_i,
182: 7.648749087422928e+00 - 4.171640244747463e+00*PETSC_i,
183: -1.829749817484586e+01 },
184: p1q2[3] = { 2.681082873627756e+00 + 3.050430199247411e+00*PETSC_i,
185: 2.681082873627756e+00 - 3.050430199247411e+00*PETSC_i,
186: 3.637834252744491e+00 },
187: p1r1[2] = { 1.000000000000000e+00 - 3.535533905932738e+00*PETSC_i,
188: 1.000000000000000e+00 + 3.535533905932738e+00*PETSC_i},
189: p1q1[2] = { 2.000000000000000e+00 + 1.414213562373095e+00*PETSC_i,
190: 2.000000000000000e+00 - 1.414213562373095e+00*PETSC_i};
191: const PetscComplex /* m == k-1 */
192: m1r5[4] = {-1.423367961376821e+02 - 1.385465094833037e+01*PETSC_i,
193: -1.423367961376821e+02 + 1.385465094833037e+01*PETSC_i,
194: 2.647367961376822e+02 - 4.814394493714596e+02*PETSC_i,
195: 2.647367961376822e+02 + 4.814394493714596e+02*PETSC_i},
196: m1q5[4] = { 5.203941240131764e+00 + 5.805856841805367e+00*PETSC_i,
197: 5.203941240131764e+00 - 5.805856841805367e+00*PETSC_i,
198: 6.796058759868242e+00 + 1.886649260140217e+00*PETSC_i,
199: 6.796058759868242e+00 - 1.886649260140217e+00*PETSC_i},
200: m1r4[3] = { 2.484269593165883e+01 + 7.460342395992306e+01*PETSC_i,
201: 2.484269593165883e+01 - 7.460342395992306e+01*PETSC_i,
202: -1.734353918633177e+02 },
203: m1q4[3] = { 4.675757014491557e+00 + 3.913489560603711e+00*PETSC_i,
204: 4.675757014491557e+00 - 3.913489560603711e+00*PETSC_i,
205: 5.648485971016893e+00 },
206: m1r3[2] = { 2.533333333333333e+01 - 2.733333333333333e+01*PETSC_i,
207: 2.533333333333333e+01 + 2.733333333333333e+01*PETSC_i},
208: m1q3[2] = { 4.000000000000000e+00 + 2.000000000000000e+00*PETSC_i,
209: 4.000000000000000e+00 - 2.000000000000000e+00*PETSC_i};
210: const PetscScalar /* m == k-1 */
211: m1remain5[2] = { 2.000000000000000e-01, 9.800000000000000e+00},
212: m1remain4[2] = {-2.500000000000000e-01, -7.750000000000000e+00},
213: m1remain3[2] = { 3.333333333333333e-01, 5.666666666666667e+00},
214: m1remain2[2] = {-0.5, -3.5},
215: remain3[4] = {1.0/6.0, 1.0/2.0, 1, 1},
216: remain2[3] = {1.0/2.0, 1, 1};
218: PetscFunctionBegin;
219: if (query) { /* query about buffer's size */
220: if (m==k+1) {
221: *remain = 0;
222: *r = *q = k+1;
223: PetscFunctionReturn(PETSC_SUCCESS); /* quick return */
224: }
225: if (m==k-1) {
226: *remain = 2;
227: if (k==5) *r = *q = 4;
228: else if (k==4) *r = *q = 3;
229: else if (k==3) *r = *q = 2;
230: else if (k==2) *r = *q = 1;
231: }
232: if (m==0) {
233: *r = *q = 0;
234: *remain = k+1;
235: }
236: } else {
237: if (m==k+1) {
238: if (k==4) {
239: for (i=0;i<5;i++) { r[i] = p1r4[i]; q[i] = p1q4[i]; }
240: } else if (k==3) {
241: for (i=0;i<4;i++) { r[i] = p1r3[i]; q[i] = p1q3[i]; }
242: } else if (k==2) {
243: for (i=0;i<3;i++) { r[i] = p1r2[i]; q[i] = p1q2[i]; }
244: } else if (k==1) {
245: for (i=0;i<2;i++) { r[i] = p1r1[i]; q[i] = p1q1[i]; }
246: }
247: PetscFunctionReturn(PETSC_SUCCESS); /* quick return */
248: }
249: if (m==k-1) {
250: if (k==5) {
251: for (i=0;i<4;i++) { r[i] = m1r5[i]; q[i] = m1q5[i]; }
252: for (i=0;i<2;i++) remain[i] = m1remain5[i];
253: } else if (k==4) {
254: for (i=0;i<3;i++) { r[i] = m1r4[i]; q[i] = m1q4[i]; }
255: for (i=0;i<2;i++) remain[i] = m1remain4[i];
256: } else if (k==3) {
257: for (i=0;i<2;i++) { r[i] = m1r3[i]; q[i] = m1q3[i]; remain[i] = m1remain3[i]; }
258: } else if (k==2) {
259: r[0] = -13.5; q[0] = 3;
260: for (i=0;i<2;i++) remain[i] = m1remain2[i];
261: }
262: }
263: if (m==0) {
264: r = q = NULL;
265: if (k==3) {
266: for (i=0;i<4;i++) remain[i] = remain3[i];
267: } else if (k==2) {
268: for (i=0;i<3;i++) remain[i] = remain2[i];
269: }
270: }
271: }
272: PetscFunctionReturn(PETSC_SUCCESS);
273: }
275: /*
276: * Product form coefficients.
277: * If query, the function returns the size necessary to store the coefficients.
278: */
279: static PetscErrorCode getcoeffsproduct(PetscInt k,PetscInt m,PetscComplex *p,PetscComplex *q,PetscComplex *mult,PetscBool query)
280: {
281: PetscInt i;
282: const PetscComplex /* m == k+1 */
283: p1p4[4] = {-5.203941240131764e+00 + 5.805856841805367e+00*PETSC_i,
284: -5.203941240131764e+00 - 5.805856841805367e+00*PETSC_i,
285: -6.796058759868242e+00 + 1.886649260140217e+00*PETSC_i,
286: -6.796058759868242e+00 - 1.886649260140217e+00*PETSC_i},
287: p1q4[5] = { 3.655694325463550e+00 + 6.543736899360086e+00*PETSC_i,
288: 3.655694325463550e+00 - 6.543736899360086e+00*PETSC_i,
289: 6.286704751729261e+00 ,
290: 5.700953298671832e+00 + 3.210265600308496e+00*PETSC_i,
291: 5.700953298671832e+00 - 3.210265600308496e+00*PETSC_i},
292: p1p3[3] = {-4.675757014491557e+00 + 3.913489560603711e+00*PETSC_i,
293: -4.675757014491557e+00 - 3.913489560603711e+00*PETSC_i,
294: -5.648485971016893e+00 },
295: p1q3[4] = { 3.212806896871536e+00 + 4.773087433276636e+00*PETSC_i,
296: 3.212806896871536e+00 - 4.773087433276636e+00*PETSC_i,
297: 4.787193103128464e+00 + 1.567476416895212e+00*PETSC_i,
298: 4.787193103128464e+00 - 1.567476416895212e+00*PETSC_i},
299: p1p2[2] = {-4.00000000000000e+00 + 2.000000000000000e+00*PETSC_i,
300: -4.00000000000000e+00 - 2.000000000000000e+00*PETSC_i},
301: p1q2[3] = { 2.681082873627756e+00 + 3.050430199247411e+00*PETSC_i,
302: 2.681082873627756e+00 - 3.050430199247411e+00*PETSC_i,
303: 3.637834252744491e+00 },
304: p1q1[2] = { 2.000000000000000e+00 + 1.414213562373095e+00*PETSC_i,
305: 2.000000000000000e+00 - 1.414213562373095e+00*PETSC_i};
306: const PetscComplex /* m == k-1 */
307: m1p5[5] = {-3.655694325463550e+00 + 6.543736899360086e+00*PETSC_i,
308: -3.655694325463550e+00 - 6.543736899360086e+00*PETSC_i,
309: -6.286704751729261e+00 ,
310: -5.700953298671832e+00 + 3.210265600308496e+00*PETSC_i,
311: -5.700953298671832e+00 - 3.210265600308496e+00*PETSC_i},
312: m1q5[4] = { 5.203941240131764e+00 + 5.805856841805367e+00*PETSC_i,
313: 5.203941240131764e+00 - 5.805856841805367e+00*PETSC_i,
314: 6.796058759868242e+00 + 1.886649260140217e+00*PETSC_i,
315: 6.796058759868242e+00 - 1.886649260140217e+00*PETSC_i},
316: m1p4[4] = {-3.212806896871536e+00 + 4.773087433276636e+00*PETSC_i,
317: -3.212806896871536e+00 - 4.773087433276636e+00*PETSC_i,
318: -4.787193103128464e+00 + 1.567476416895212e+00*PETSC_i,
319: -4.787193103128464e+00 - 1.567476416895212e+00*PETSC_i},
320: m1q4[3] = { 4.675757014491557e+00 + 3.913489560603711e+00*PETSC_i,
321: 4.675757014491557e+00 - 3.913489560603711e+00*PETSC_i,
322: 5.648485971016893e+00 },
323: m1p3[3] = {-2.681082873627756e+00 + 3.050430199247411e+00*PETSC_i,
324: -2.681082873627756e+00 - 3.050430199247411e+00*PETSC_i,
325: -3.637834252744491e+00 },
326: m1q3[2] = { 4.000000000000000e+00 + 2.000000000000000e+00*PETSC_i,
327: 4.000000000000000e+00 - 2.000000000000001e+00*PETSC_i},
328: m1p2[2] = {-2.000000000000000e+00 + 1.414213562373095e+00*PETSC_i,
329: -2.000000000000000e+00 - 1.414213562373095e+00*PETSC_i};
331: PetscFunctionBegin;
332: if (query) {
333: if (m == k+1) {
334: *mult = 1;
335: *p = k;
336: *q = k+1;
337: PetscFunctionReturn(PETSC_SUCCESS);
338: }
339: if (m==k-1) {
340: *mult = 1;
341: *p = k;
342: *q = k-1;
343: }
344: } else {
345: if (m == k+1) {
346: *mult = PetscPowInt(-1,m);
347: *mult *= m;
348: if (k==4) {
349: for (i=0;i<4;i++) { p[i] = p1p4[i]; q[i] = p1q4[i]; }
350: q[4] = p1q4[4];
351: } else if (k==3) {
352: for (i=0;i<3;i++) { p[i] = p1p3[i]; q[i] = p1q3[i]; }
353: q[3] = p1q3[3];
354: } else if (k==2) {
355: for (i=0;i<2;i++) { p[i] = p1p2[i]; q[i] = p1q2[i]; }
356: q[2] = p1q2[2];
357: } else if (k==1) {
358: p[0] = -3;
359: for (i=0;i<2;i++) q[i] = p1q1[i];
360: }
361: PetscFunctionReturn(PETSC_SUCCESS);
362: }
363: if (m==k-1) {
364: *mult = PetscPowInt(-1,m);
365: *mult /= k;
366: if (k==5) {
367: for (i=0;i<4;i++) { p[i] = m1p5[i]; q[i] = m1q5[i]; }
368: p[4] = m1p5[4];
369: } else if (k==4) {
370: for (i=0;i<3;i++) { p[i] = m1p4[i]; q[i] = m1q4[i]; }
371: p[3] = m1p4[3];
372: } else if (k==3) {
373: for (i=0;i<2;i++) { p[i] = m1p3[i]; q[i] = m1q3[i]; }
374: p[2] = m1p3[2];
375: } else if (k==2) {
376: for (i=0;i<2;i++) p[i] = m1p2[i];
377: q[0] = 3;
378: }
379: }
380: }
381: PetscFunctionReturn(PETSC_SUCCESS);
382: }
383: #endif /* PETSC_HAVE_COMPLEX */
385: #if PetscDefined(USE_COMPLEX)
386: static PetscErrorCode getisreal(PetscInt n,PetscComplex *a,PetscBool *result)
387: {
388: PetscInt i;
390: PetscFunctionBegin;
391: *result=PETSC_TRUE;
392: for (i=0;i<n&&*result;i++) {
393: if (PetscImaginaryPartComplex(a[i])) *result=PETSC_FALSE;
394: }
395: PetscFunctionReturn(PETSC_SUCCESS);
396: }
397: #endif
399: /*
400: * Matrix exponential implementation based on algorithm and matlab code by Stefan Guettel
401: * and Yuji Nakatsukasa
402: *
403: * Stefan Guettel and Yuji Nakatsukasa, "Scaled and Squared Subdiagonal Pade
404: * Approximation for the Matrix Exponential",
405: * SIAM J. Matrix Anal. Appl. 37(1):145-170, 2016.
406: * https://doi.org/10.1137/15M1027553
407: */
408: static PetscErrorCode FNEvaluateFunctionMat_Exp_GuettelNakatsukasa(FN fn,Mat A,Mat B)
409: {
410: #if !PetscDefined(HAVE_COMPLEX)
411: PetscFunctionBegin;
412: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"This function requires C99 or C++ complex support");
413: #else
414: PetscInt i,j,n_,s,k,m,mod;
415: PetscBLASInt n=0,n2=0,irsize=0,rsizediv2,ipsize=0,iremainsize=0,*piv,minlen,lwork=0,one=1;
416: PetscReal nrm,shift=0.0;
417: #if PetscDefined(USE_COMPLEX)
418: PetscReal *rwork=NULL;
419: #endif
420: PetscComplex *As,*RR,*RR2,*expmA,*expmA2,*Maux,*Maux2,rsize,*r,psize,*p,remainsize,*remainterm,*rootp,*rootq,mult=0.0,scale,cone=1.0,czero=0.0,*aux;
421: PetscScalar *Ba,*Ba2,*sMaux,*wr,*wi,expshift,sone=1.0,szero=0.0,*saux;
422: const PetscScalar *Aa;
423: PetscBool isreal,flg;
424: PetscBLASInt query=-1;
425: PetscScalar work1,*work;
427: PetscFunctionBegin;
428: PetscCall(MatGetSize(A,&n_,NULL));
429: PetscCall(PetscBLASIntCast(n_,&n));
430: PetscCall(MatDenseGetArrayRead(A,&Aa));
431: PetscCall(MatDenseGetArray(B,&Ba));
432: Ba2 = Ba;
433: PetscCall(PetscBLASIntCast(n*n,&n2));
435: PetscCall(PetscMalloc2(n2,&sMaux,n2,&Maux));
436: Maux2 = Maux;
437: PetscCall(PetscOptionsGetReal(NULL,NULL,"-fn_expm_estimated_eig",&shift,&flg));
438: if (!flg) {
439: PetscCall(PetscMalloc2(n,&wr,n,&wi));
440: PetscCall(PetscArraycpy(sMaux,Aa,n2));
441: /* estimate rightmost eigenvalue and shift A with it */
442: #if !PetscDefined(USE_COMPLEX)
443: PetscCallLAPACKInfo("LAPACKgeev",LAPACKgeev_("N","N",&n,sMaux,&n,wr,wi,NULL,&n,NULL,&n,&work1,&query,&info));
444: PetscCall(PetscBLASIntCast((PetscInt)work1,&lwork));
445: PetscCall(PetscMalloc1(lwork,&work));
446: PetscCallLAPACKInfo("LAPACKgeev",LAPACKgeev_("N","N",&n,sMaux,&n,wr,wi,NULL,&n,NULL,&n,work,&lwork,&info));
447: PetscCall(PetscFree(work));
448: #else
449: PetscCall(PetscArraycpy(Maux,Aa,n2));
450: PetscCallLAPACKInfo("LAPACKgeev",LAPACKgeev_("N","N",&n,Maux,&n,wr,NULL,&n,NULL,&n,&work1,&query,rwork,&info));
451: PetscCall(PetscBLASIntCast((PetscInt)PetscRealPart(work1),&lwork));
452: PetscCall(PetscMalloc2(2*n,&rwork,lwork,&work));
453: PetscCallLAPACKInfo("LAPACKgeev",LAPACKgeev_("N","N",&n,Maux,&n,wr,NULL,&n,NULL,&n,work,&lwork,rwork,&info));
454: PetscCall(PetscFree2(rwork,work));
455: #endif
456: PetscCall(PetscLogFlops(25.0*n*n*n+(n*n*n)/3.0+1.0*n*n*n));
458: shift = PetscRealPart(wr[0]);
459: for (i=1;i<n;i++) {
460: if (PetscRealPart(wr[i]) > shift) shift = PetscRealPart(wr[i]);
461: }
462: PetscCall(PetscFree2(wr,wi));
463: }
464: /* shift so that largest real part is (about) 0 */
465: PetscCall(PetscArraycpy(sMaux,Aa,n2));
466: if (shift) {
467: for (i=0;i<n;i++) sMaux[i+i*n] -= shift;
468: PetscCall(PetscLogFlops(1.0*n));
469: }
470: #if PetscDefined(USE_COMPLEX)
471: PetscCall(PetscArraycpy(Maux,Aa,n2));
472: if (shift) {
473: for (i=0;i<n;i++) Maux[i+i*n] -= shift;
474: PetscCall(PetscLogFlops(1.0*n));
475: }
476: #endif
478: /* estimate norm(A) and select the scaling factor */
479: nrm = LAPACKlange_("O",&n,&n,sMaux,&n,NULL);
480: PetscCall(PetscLogFlops(1.0*n*n));
481: PetscCall(sexpm_params(nrm,&s,&k,&m));
482: if (s==0 && k==1 && m==0) { /* exp(A) = I+A to eps! */
483: if (shift) expshift = PetscExpReal(shift);
484: for (i=0;i<n;i++) sMaux[i+i*n] += 1.0;
485: if (shift) {
486: PetscCallBLAS("BLASscal",BLASscal_(&n2,&expshift,sMaux,&one));
487: PetscCall(PetscLogFlops(1.0*(n+n2)));
488: } else PetscCall(PetscLogFlops(1.0*n));
489: PetscCall(PetscArraycpy(Ba,sMaux,n2));
490: PetscCall(PetscFree2(sMaux,Maux));
491: PetscCall(MatDenseRestoreArrayRead(A,&Aa));
492: PetscCall(MatDenseRestoreArray(B,&Ba));
493: PetscFunctionReturn(PETSC_SUCCESS); /* quick return */
494: }
496: PetscCall(PetscMalloc4(n2,&expmA,n2,&As,n2,&RR,n,&piv));
497: expmA2 = expmA; RR2 = RR;
498: /* scale matrix */
499: #if !PetscDefined(USE_COMPLEX)
500: for (i=0;i<n2;i++) {
501: As[i] = sMaux[i];
502: }
503: #else
504: PetscCall(PetscArraycpy(As,sMaux,n2));
505: #endif
506: scale = 1.0/PetscPowRealInt(2.0,s);
507: PetscCallBLAS("BLASCOMPLEXscal",BLASCOMPLEXscal_(&n2,&scale,As,&one));
508: PetscCall(SlepcLogFlopsComplex(1.0*n2));
510: /* evaluate Pade approximant (partial fraction or product form) */
511: if (fn->method==3 || !m) { /* partial fraction */
512: PetscCall(getcoeffs(k,m,&rsize,&psize,&remainsize,PETSC_TRUE));
513: PetscCall(PetscBLASIntCast((PetscInt)PetscRealPartComplex(rsize),&irsize));
514: PetscCall(PetscBLASIntCast((PetscInt)PetscRealPartComplex(psize),&ipsize));
515: PetscCall(PetscBLASIntCast((PetscInt)PetscRealPartComplex(remainsize),&iremainsize));
516: PetscCall(PetscMalloc3(irsize,&r,ipsize,&p,iremainsize,&remainterm));
517: PetscCall(getcoeffs(k,m,r,p,remainterm,PETSC_FALSE));
519: PetscCall(PetscArrayzero(expmA,n2));
520: #if !PetscDefined(USE_COMPLEX)
521: isreal = PETSC_TRUE;
522: #else
523: PetscCall(getisreal(n2,Maux,&isreal));
524: #endif
525: if (isreal) {
526: rsizediv2 = irsize/2;
527: for (i=0;i<rsizediv2;i++) { /* use partial fraction to get R(As) */
528: PetscCall(PetscArraycpy(Maux,As,n2));
529: PetscCall(PetscArrayzero(RR,n2));
530: for (j=0;j<n;j++) {
531: Maux[j+j*n] -= p[2*i];
532: RR[j+j*n] = r[2*i];
533: }
534: PetscCallLAPACKInfo("LAPACKCOMPLEXgesv",LAPACKCOMPLEXgesv_(&n,&n,Maux,&n,piv,RR,&n,&info));
535: for (j=0;j<n2;j++) {
536: expmA[j] += RR[j] + PetscConj(RR[j]);
537: }
538: /* loop(n) + gesv + loop(n2) */
539: PetscCall(SlepcLogFlopsComplex(1.0*n+(2.0*n*n*n/3.0+2.0*n*n*n)+2.0*n2));
540: }
542: mod = ipsize % 2;
543: if (mod) {
544: PetscCall(PetscArraycpy(Maux,As,n2));
545: PetscCall(PetscArrayzero(RR,n2));
546: for (j=0;j<n;j++) {
547: Maux[j+j*n] -= p[ipsize-1];
548: RR[j+j*n] = r[irsize-1];
549: }
550: PetscCallLAPACKInfo("LAPACKCOMPLEXgesv",LAPACKCOMPLEXgesv_(&n,&n,Maux,&n,piv,RR,&n,&info));
551: for (j=0;j<n2;j++) {
552: expmA[j] += RR[j];
553: }
554: PetscCall(SlepcLogFlopsComplex(1.0*n+(2.0*n*n*n/3.0+2.0*n*n*n)+1.0*n2));
555: }
556: } else { /* complex */
557: for (i=0;i<irsize;i++) { /* use partial fraction to get R(As) */
558: PetscCall(PetscArraycpy(Maux,As,n2));
559: PetscCall(PetscArrayzero(RR,n2));
560: for (j=0;j<n;j++) {
561: Maux[j+j*n] -= p[i];
562: RR[j+j*n] = r[i];
563: }
564: PetscCallLAPACKInfo("LAPACKCOMPLEXgesv",LAPACKCOMPLEXgesv_(&n,&n,Maux,&n,piv,RR,&n,&info));
565: for (j=0;j<n2;j++) {
566: expmA[j] += RR[j];
567: }
568: PetscCall(SlepcLogFlopsComplex(1.0*n+(2.0*n*n*n/3.0+2.0*n*n*n)+1.0*n2));
569: }
570: }
571: for (i=0;i<iremainsize;i++) {
572: if (!i) {
573: PetscCall(PetscArrayzero(RR,n2));
574: for (j=0;j<n;j++) {
575: RR[j+j*n] = remainterm[iremainsize-1];
576: }
577: } else {
578: PetscCall(PetscArraycpy(RR,As,n2));
579: for (j=1;j<i;j++) {
580: PetscCallBLAS("BLASCOMPLEXgemm",BLASCOMPLEXgemm_("N","N",&n,&n,&n,&cone,RR,&n,RR,&n,&czero,Maux,&n));
581: SlepcSwap(RR,Maux,aux);
582: PetscCall(SlepcLogFlopsComplex(2.0*n*n*n));
583: }
584: PetscCallBLAS("BLASCOMPLEXscal",BLASCOMPLEXscal_(&n2,&remainterm[iremainsize-1-i],RR,&one));
585: PetscCall(SlepcLogFlopsComplex(1.0*n2));
586: }
587: for (j=0;j<n2;j++) {
588: expmA[j] += RR[j];
589: }
590: PetscCall(SlepcLogFlopsComplex(1.0*n2));
591: }
592: PetscCall(PetscFree3(r,p,remainterm));
593: } else { /* product form, default */
594: PetscCall(getcoeffsproduct(k,m,&rsize,&psize,&mult,PETSC_TRUE));
595: PetscCall(PetscBLASIntCast((PetscInt)PetscRealPartComplex(rsize),&irsize));
596: PetscCall(PetscBLASIntCast((PetscInt)PetscRealPartComplex(psize),&ipsize));
597: PetscCall(PetscMalloc2(irsize,&rootp,ipsize,&rootq));
598: PetscCall(getcoeffsproduct(k,m,rootp,rootq,&mult,PETSC_FALSE));
600: PetscCall(PetscArrayzero(expmA,n2));
601: for (i=0;i<n;i++) { /* initialize */
602: expmA[i+i*n] = 1.0;
603: }
604: minlen = PetscMin(irsize,ipsize);
605: for (i=0;i<minlen;i++) {
606: PetscCall(PetscArraycpy(RR,As,n2));
607: for (j=0;j<n;j++) {
608: RR[j+j*n] -= rootp[i];
609: }
610: PetscCallBLAS("BLASCOMPLEXgemm",BLASCOMPLEXgemm_("N","N",&n,&n,&n,&cone,RR,&n,expmA,&n,&czero,Maux,&n));
611: SlepcSwap(expmA,Maux,aux);
612: PetscCall(PetscArraycpy(RR,As,n2));
613: for (j=0;j<n;j++) {
614: RR[j+j*n] -= rootq[i];
615: }
616: PetscCallLAPACKInfo("LAPACKCOMPLEXgesv",LAPACKCOMPLEXgesv_(&n,&n,RR,&n,piv,expmA,&n,&info));
617: /* loop(n) + gemm + loop(n) + gesv */
618: PetscCall(SlepcLogFlopsComplex(1.0*n+(2.0*n*n*n)+1.0*n+(2.0*n*n*n/3.0+2.0*n*n*n)));
619: }
620: /* extra numerator */
621: for (i=minlen;i<irsize;i++) {
622: PetscCall(PetscArraycpy(RR,As,n2));
623: for (j=0;j<n;j++) {
624: RR[j+j*n] -= rootp[i];
625: }
626: PetscCallBLAS("BLASCOMPLEXgemm",BLASCOMPLEXgemm_("N","N",&n,&n,&n,&cone,RR,&n,expmA,&n,&czero,Maux,&n));
627: SlepcSwap(expmA,Maux,aux);
628: PetscCall(SlepcLogFlopsComplex(1.0*n+2.0*n*n*n));
629: }
630: /* extra denominator */
631: for (i=minlen;i<ipsize;i++) {
632: PetscCall(PetscArraycpy(RR,As,n2));
633: for (j=0;j<n;j++) RR[j+j*n] -= rootq[i];
634: PetscCallLAPACKInfo("LAPACKCOMPLEXgesv",LAPACKCOMPLEXgesv_(&n,&n,RR,&n,piv,expmA,&n,&info));
635: PetscCall(SlepcLogFlopsComplex(1.0*n+(2.0*n*n*n/3.0+2.0*n*n*n)));
636: }
637: PetscCallBLAS("BLASCOMPLEXscal",BLASCOMPLEXscal_(&n2,&mult,expmA,&one));
638: PetscCall(SlepcLogFlopsComplex(1.0*n2));
639: PetscCall(PetscFree2(rootp,rootq));
640: }
642: #if !PetscDefined(USE_COMPLEX)
643: for (i=0;i<n2;i++) {
644: Ba2[i] = PetscRealPartComplex(expmA[i]);
645: }
646: #else
647: PetscCall(PetscArraycpy(Ba2,expmA,n2));
648: #endif
650: /* perform repeated squaring */
651: for (i=0;i<s;i++) { /* final squaring */
652: PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&n,&n,&n,&sone,Ba2,&n,Ba2,&n,&szero,sMaux,&n));
653: SlepcSwap(Ba2,sMaux,saux);
654: PetscCall(PetscLogFlops(2.0*n*n*n));
655: }
656: if (Ba2!=Ba) {
657: PetscCall(PetscArraycpy(Ba,Ba2,n2));
658: sMaux = Ba2;
659: }
660: if (shift) {
661: expshift = PetscExpReal(shift);
662: PetscCallBLAS("BLASscal",BLASscal_(&n2,&expshift,Ba,&one));
663: PetscCall(PetscLogFlops(1.0*n2));
664: }
666: /* restore pointers */
667: Maux = Maux2; expmA = expmA2; RR = RR2;
668: PetscCall(PetscFree2(sMaux,Maux));
669: PetscCall(PetscFree4(expmA,As,RR,piv));
670: PetscCall(MatDenseRestoreArrayRead(A,&Aa));
671: PetscCall(MatDenseRestoreArray(B,&Ba));
672: PetscFunctionReturn(PETSC_SUCCESS);
673: #endif
674: }
676: #define SMALLN 100
678: /*
679: * Function needed to compute optimal parameters (required workspace is 3*n*n)
680: */
681: static PetscInt ell(PetscBLASInt n,PetscScalar *A,PetscReal coeff,PetscInt m,PetscScalar *work,PetscRandom rand)
682: {
683: PetscScalar *Ascaled=work;
684: PetscReal nrm,alpha,beta,rwork[1];
685: PetscInt t;
686: PetscBLASInt i,j;
688: PetscFunctionBegin;
689: beta = PetscPowReal(coeff,1.0/(2*m+1));
690: for (i=0;i<n;i++)
691: for (j=0;j<n;j++)
692: Ascaled[i+j*n] = beta*PetscAbsScalar(A[i+j*n]);
693: nrm = LAPACKlange_("O",&n,&n,A,&n,rwork);
694: PetscCall(PetscLogFlops(2.0*n*n));
695: PetscCall(SlepcNormAm(n,Ascaled,2*m+1,work+n*n,rand,&alpha));
696: alpha /= nrm;
697: t = PetscMax((PetscInt)PetscCeilReal(PetscLogReal(2.0*alpha/PETSC_MACHINE_EPSILON)/PetscLogReal(2.0)/(2*m)),0);
698: PetscFunctionReturn(t);
699: }
701: /*
702: * Compute scaling parameter (s) and order of Pade approximant (m) (required workspace is 4*n*n)
703: */
704: static PetscErrorCode expm_params(PetscInt n,PetscScalar **Apowers,PetscInt *s,PetscInt *m,PetscScalar *work)
705: {
706: PetscScalar sfactor,sone=1.0,szero=0.0,*A=Apowers[0],*Ascaled;
707: PetscReal d4,d6,d8,d10,eta1,eta3,eta4,eta5,rwork[1];
708: PetscBLASInt n_=0,n2,one=1;
709: PetscRandom rand;
710: const PetscReal coeff[5] = { 9.92063492063492e-06, 9.94131285136576e-11, /* backward error function */
711: 2.22819456055356e-16, 1.69079293431187e-22, 8.82996160201868e-36 };
712: const PetscReal theta[5] = { 1.495585217958292e-002, /* m = 3 */
713: 2.539398330063230e-001, /* m = 5 */
714: 9.504178996162932e-001, /* m = 7 */
715: 2.097847961257068e+000, /* m = 9 */
716: 5.371920351148152e+000 }; /* m = 13 */
718: PetscFunctionBegin;
719: *s = 0;
720: *m = 13;
721: PetscCall(PetscBLASIntCast(n,&n_));
722: PetscCall(PetscRandomCreate(PETSC_COMM_SELF,&rand));
723: d4 = PetscPowReal(LAPACKlange_("O",&n_,&n_,Apowers[2],&n_,rwork),1.0/4.0);
724: if (d4==0.0) { /* safeguard for the case A = 0 */
725: *m = 3;
726: goto done;
727: }
728: d6 = PetscPowReal(LAPACKlange_("O",&n_,&n_,Apowers[3],&n_,rwork),1.0/6.0);
729: PetscCall(PetscLogFlops(2.0*n*n));
730: eta1 = PetscMax(d4,d6);
731: if (eta1<=theta[0] && !ell(n_,A,coeff[0],3,work,rand)) {
732: *m = 3;
733: goto done;
734: }
735: if (eta1<=theta[1] && !ell(n_,A,coeff[1],5,work,rand)) {
736: *m = 5;
737: goto done;
738: }
739: if (n<SMALLN) {
740: PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&n_,&n_,&n_,&sone,Apowers[2],&n_,Apowers[2],&n_,&szero,work,&n_));
741: d8 = PetscPowReal(LAPACKlange_("O",&n_,&n_,work,&n_,rwork),1.0/8.0);
742: PetscCall(PetscLogFlops(2.0*n*n*n+1.0*n*n));
743: } else {
744: PetscCall(SlepcNormAm(n_,Apowers[2],2,work,rand,&d8));
745: d8 = PetscPowReal(d8,1.0/8.0);
746: }
747: eta3 = PetscMax(d6,d8);
748: if (eta3<=theta[2] && !ell(n_,A,coeff[2],7,work,rand)) {
749: *m = 7;
750: goto done;
751: }
752: if (eta3<=theta[3] && !ell(n_,A,coeff[3],9,work,rand)) {
753: *m = 9;
754: goto done;
755: }
756: if (n<SMALLN) {
757: PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&n_,&n_,&n_,&sone,Apowers[2],&n_,Apowers[3],&n_,&szero,work,&n_));
758: d10 = PetscPowReal(LAPACKlange_("O",&n_,&n_,work,&n_,rwork),1.0/10.0);
759: PetscCall(PetscLogFlops(2.0*n*n*n+1.0*n*n));
760: } else {
761: PetscCall(SlepcNormAm(n_,Apowers[1],5,work,rand,&d10));
762: d10 = PetscPowReal(d10,1.0/10.0);
763: }
764: eta4 = PetscMax(d8,d10);
765: eta5 = PetscMin(eta3,eta4);
766: *s = PetscMax((PetscInt)PetscCeilReal(PetscLogReal(eta5/theta[4])/PetscLogReal(2.0)),0);
767: if (*s) {
768: Ascaled = work+3*n*n;
769: n2 = n_*n_;
770: PetscCallBLAS("BLAScopy",BLAScopy_(&n2,A,&one,Ascaled,&one));
771: sfactor = PetscPowRealInt(2.0,-(*s));
772: PetscCallBLAS("BLASscal",BLASscal_(&n2,&sfactor,Ascaled,&one));
773: PetscCall(PetscLogFlops(1.0*n*n));
774: } else Ascaled = A;
775: *s += ell(n_,Ascaled,coeff[4],13,work,rand);
776: done:
777: PetscCall(PetscRandomDestroy(&rand));
778: PetscFunctionReturn(PETSC_SUCCESS);
779: }
781: /*
782: * Matrix exponential implementation based on algorithm and matlab code by N. Higham and co-authors
783: *
784: * N. J. Higham, "The scaling and squaring method for the matrix exponential
785: * revisited", SIAM J. Matrix Anal. Appl. 26(4):1179-1193, 2005.
786: */
787: PetscErrorCode FNEvaluateFunctionMat_Exp_Higham(FN fn,Mat A,Mat B)
788: {
789: PetscBLASInt n_=0,n2,*ipiv,one=1;
790: PetscInt n,m,j,s;
791: PetscScalar scale,smone=-1.0,sone=1.0,stwo=2.0,szero=0.0;
792: PetscScalar *Ba,*Apowers[5],*Q,*P,*W,*work,*aux;
793: const PetscScalar *Aa,*c;
794: const PetscScalar c3[4] = { 120, 60, 12, 1 };
795: const PetscScalar c5[6] = { 30240, 15120, 3360, 420, 30, 1 };
796: const PetscScalar c7[8] = { 17297280, 8648640, 1995840, 277200, 25200, 1512, 56, 1 };
797: const PetscScalar c9[10] = { 17643225600.0, 8821612800.0, 2075673600, 302702400, 30270240,
798: 2162160, 110880, 3960, 90, 1 };
799: const PetscScalar c13[14] = { 64764752532480000.0, 32382376266240000.0, 7771770303897600.0,
800: 1187353796428800.0, 129060195264000.0, 10559470521600.0,
801: 670442572800.0, 33522128640.0, 1323241920.0,
802: 40840800, 960960, 16380, 182, 1 };
804: PetscFunctionBegin;
805: PetscCall(MatDenseGetArrayRead(A,&Aa));
806: PetscCall(MatDenseGetArray(B,&Ba));
807: PetscCall(MatGetSize(A,&n,NULL));
808: PetscCall(PetscBLASIntCast(n,&n_));
809: n2 = n_*n_;
810: PetscCall(PetscMalloc2(8*n*n,&work,n,&ipiv));
812: /* Matrix powers */
813: Apowers[0] = work; /* Apowers[0] = A */
814: Apowers[1] = Apowers[0] + n*n; /* Apowers[1] = A^2 */
815: Apowers[2] = Apowers[1] + n*n; /* Apowers[2] = A^4 */
816: Apowers[3] = Apowers[2] + n*n; /* Apowers[3] = A^6 */
817: Apowers[4] = Apowers[3] + n*n; /* Apowers[4] = A^8 */
819: PetscCall(PetscArraycpy(Apowers[0],Aa,n2));
820: PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&n_,&n_,&n_,&sone,Apowers[0],&n_,Apowers[0],&n_,&szero,Apowers[1],&n_));
821: PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&n_,&n_,&n_,&sone,Apowers[1],&n_,Apowers[1],&n_,&szero,Apowers[2],&n_));
822: PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&n_,&n_,&n_,&sone,Apowers[1],&n_,Apowers[2],&n_,&szero,Apowers[3],&n_));
823: PetscCall(PetscLogFlops(6.0*n*n*n));
825: /* Compute scaling parameter and order of Pade approximant */
826: PetscCall(expm_params(n,Apowers,&s,&m,Apowers[4]));
828: if (s) { /* rescale */
829: for (j=0;j<4;j++) {
830: scale = PetscPowRealInt(2.0,-PetscMax(2*j,1)*s);
831: PetscCallBLAS("BLASscal",BLASscal_(&n2,&scale,Apowers[j],&one));
832: }
833: PetscCall(PetscLogFlops(4.0*n*n));
834: }
836: /* Evaluate the Pade approximant */
837: switch (m) {
838: case 3: c = c3; break;
839: case 5: c = c5; break;
840: case 7: c = c7; break;
841: case 9: c = c9; break;
842: case 13: c = c13; break;
843: default: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Wrong value of m %" PetscInt_FMT,m);
844: }
845: P = Ba;
846: Q = Apowers[4] + n*n;
847: W = Q + n*n;
848: switch (m) {
849: case 3:
850: case 5:
851: case 7:
852: case 9:
853: if (m==9) PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&n_,&n_,&n_,&sone,Apowers[1],&n_,Apowers[3],&n_,&szero,Apowers[4],&n_));
854: PetscCall(PetscArrayzero(P,n2));
855: PetscCall(PetscArrayzero(Q,n2));
856: for (j=0;j<n;j++) {
857: P[j+j*n] = c[1];
858: Q[j+j*n] = c[0];
859: }
860: for (j=m;j>=3;j-=2) {
861: PetscCallBLAS("BLASaxpy",BLASaxpy_(&n2,&c[j],Apowers[(j+1)/2-1],&one,P,&one));
862: PetscCallBLAS("BLASaxpy",BLASaxpy_(&n2,&c[j-1],Apowers[(j+1)/2-1],&one,Q,&one));
863: PetscCall(PetscLogFlops(4.0*n*n));
864: }
865: PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&n_,&n_,&n_,&sone,Apowers[0],&n_,P,&n_,&szero,W,&n_));
866: PetscCall(PetscLogFlops(2.0*n*n*n));
867: SlepcSwap(P,W,aux);
868: break;
869: case 13:
870: /* P = A*(Apowers[3]*(c[13]*Apowers[3] + c[11]*Apowers[2] + c[9]*Apowers[1])
871: + c[7]*Apowers[3] + c[5]*Apowers[2] + c[3]*Apowers[1] + c[1]*I) */
872: PetscCallBLAS("BLAScopy",BLAScopy_(&n2,Apowers[3],&one,P,&one));
873: PetscCallBLAS("BLASscal",BLASscal_(&n2,&c[13],P,&one));
874: PetscCallBLAS("BLASaxpy",BLASaxpy_(&n2,&c[11],Apowers[2],&one,P,&one));
875: PetscCallBLAS("BLASaxpy",BLASaxpy_(&n2,&c[9],Apowers[1],&one,P,&one));
876: PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&n_,&n_,&n_,&sone,Apowers[3],&n_,P,&n_,&szero,W,&n_));
877: PetscCall(PetscLogFlops(5.0*n*n+2.0*n*n*n));
878: PetscCall(PetscArrayzero(P,n2));
879: for (j=0;j<n;j++) P[j+j*n] = c[1];
880: PetscCallBLAS("BLASaxpy",BLASaxpy_(&n2,&c[7],Apowers[3],&one,P,&one));
881: PetscCallBLAS("BLASaxpy",BLASaxpy_(&n2,&c[5],Apowers[2],&one,P,&one));
882: PetscCallBLAS("BLASaxpy",BLASaxpy_(&n2,&c[3],Apowers[1],&one,P,&one));
883: PetscCallBLAS("BLASaxpy",BLASaxpy_(&n2,&sone,P,&one,W,&one));
884: PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&n_,&n_,&n_,&sone,Apowers[0],&n_,W,&n_,&szero,P,&n_));
885: PetscCall(PetscLogFlops(7.0*n*n+2.0*n*n*n));
886: /* Q = Apowers[3]*(c[12]*Apowers[3] + c[10]*Apowers[2] + c[8]*Apowers[1])
887: + c[6]*Apowers[3] + c[4]*Apowers[2] + c[2]*Apowers[1] + c[0]*I */
888: PetscCallBLAS("BLAScopy",BLAScopy_(&n2,Apowers[3],&one,Q,&one));
889: PetscCallBLAS("BLASscal",BLASscal_(&n2,&c[12],Q,&one));
890: PetscCallBLAS("BLASaxpy",BLASaxpy_(&n2,&c[10],Apowers[2],&one,Q,&one));
891: PetscCallBLAS("BLASaxpy",BLASaxpy_(&n2,&c[8],Apowers[1],&one,Q,&one));
892: PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&n_,&n_,&n_,&sone,Apowers[3],&n_,Q,&n_,&szero,W,&n_));
893: PetscCall(PetscLogFlops(5.0*n*n+2.0*n*n*n));
894: PetscCall(PetscArrayzero(Q,n2));
895: for (j=0;j<n;j++) Q[j+j*n] = c[0];
896: PetscCallBLAS("BLASaxpy",BLASaxpy_(&n2,&c[6],Apowers[3],&one,Q,&one));
897: PetscCallBLAS("BLASaxpy",BLASaxpy_(&n2,&c[4],Apowers[2],&one,Q,&one));
898: PetscCallBLAS("BLASaxpy",BLASaxpy_(&n2,&c[2],Apowers[1],&one,Q,&one));
899: PetscCallBLAS("BLASaxpy",BLASaxpy_(&n2,&sone,W,&one,Q,&one));
900: PetscCall(PetscLogFlops(7.0*n*n));
901: break;
902: default: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Wrong value of m %" PetscInt_FMT,m);
903: }
904: PetscCallBLAS("BLASaxpy",BLASaxpy_(&n2,&smone,P,&one,Q,&one));
905: PetscCallLAPACKInfo("LAPACKgesv",LAPACKgesv_(&n_,&n_,Q,&n_,ipiv,P,&n_,&info));
906: PetscCallBLAS("BLASscal",BLASscal_(&n2,&stwo,P,&one));
907: for (j=0;j<n;j++) P[j+j*n] += 1.0;
908: PetscCall(PetscLogFlops(2.0*n*n*n/3.0+4.0*n*n));
910: /* Squaring */
911: for (j=1;j<=s;j++) {
912: PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&n_,&n_,&n_,&sone,P,&n_,P,&n_,&szero,W,&n_));
913: SlepcSwap(P,W,aux);
914: }
915: if (P!=Ba) PetscCall(PetscArraycpy(Ba,P,n2));
916: PetscCall(PetscLogFlops(2.0*n*n*n*s));
918: PetscCall(PetscFree2(work,ipiv));
919: PetscCall(MatDenseRestoreArrayRead(A,&Aa));
920: PetscCall(MatDenseRestoreArray(B,&Ba));
921: PetscFunctionReturn(PETSC_SUCCESS);
922: }
924: #if PetscDefined(HAVE_CUDA)
925: #include "../src/sys/classes/fn/impls/cuda/fnutilcuda.h"
926: #include <slepccupmblas.h>
928: PetscErrorCode FNEvaluateFunctionMat_Exp_Pade_CUDA(FN fn,Mat A,Mat B)
929: {
930: PetscBLASInt n=0,ld,ld2,*d_ipiv,*d_info,info,one=1;
931: PetscInt m,k,sexp;
932: PetscBool odd;
933: const PetscInt p=MAX_PADE;
934: PetscReal c[MAX_PADE+1],s;
935: PetscScalar scale,smone=-1.0,sone=1.0,stwo=2.0,szero=0.0;
936: const PetscScalar *Aa;
937: PetscScalar *d_Ba,*d_As,*d_A2,*d_Q,*d_P,*d_W,*aux,**ppP,**d_ppP,**ppQ,**d_ppQ;
938: cublasHandle_t cublasv2handle;
940: PetscFunctionBegin;
941: PetscCall(PetscDeviceInitialize(PETSC_DEVICE_CUDA)); /* For CUDA event timers */
942: PetscCall(PetscCUBLASGetHandle(&cublasv2handle));
943: PetscCall(MatGetSize(A,&m,NULL));
944: PetscCall(PetscBLASIntCast(m,&n));
945: ld = n;
946: ld2 = ld*ld;
947: if (A==B) {
948: PetscCallCUDA(cudaMalloc((void **)&d_As,sizeof(PetscScalar)*m*m));
949: PetscCall(MatDenseCUDAGetArrayRead(A,&Aa));
950: PetscCallCUDA(cudaMemcpy(d_As,Aa,sizeof(PetscScalar)*ld2,cudaMemcpyDeviceToDevice));
951: PetscCall(MatDenseCUDARestoreArrayRead(A,&Aa));
952: } else PetscCall(MatDenseCUDAGetArrayRead(A,(const PetscScalar**)&d_As));
953: PetscCall(MatDenseCUDAGetArrayWrite(B,&d_Ba));
955: PetscCallCUDA(cudaMalloc((void **)&d_Q,sizeof(PetscScalar)*m*m));
956: PetscCallCUDA(cudaMalloc((void **)&d_W,sizeof(PetscScalar)*m*m));
957: PetscCallCUDA(cudaMalloc((void **)&d_A2,sizeof(PetscScalar)*m*m));
958: PetscCallCUDA(cudaMalloc((void **)&d_ipiv,sizeof(PetscBLASInt)*ld));
959: PetscCallCUDA(cudaMalloc((void **)&d_info,sizeof(PetscBLASInt)));
960: PetscCallCUDA(cudaMalloc((void **)&d_ppP,sizeof(PetscScalar*)));
961: PetscCallCUDA(cudaMalloc((void **)&d_ppQ,sizeof(PetscScalar*)));
963: PetscCall(PetscMalloc1(1,&ppP));
964: PetscCall(PetscMalloc1(1,&ppQ));
966: d_P = d_Ba;
967: PetscCall(PetscLogGpuTimeBegin());
969: /* Pade' coefficients */
970: c[0] = 1.0;
971: for (k=1;k<=p;k++) c[k] = c[k-1]*(p+1-k)/(k*(2*p+1-k));
973: /* Scaling */
974: PetscCallCUBLAS(cublasXnrm2(cublasv2handle,ld2,d_As,one,&s));
975: if (s>0.5) {
976: sexp = PetscMax(0,(int)(PetscLogReal(s)/PetscLogReal(2.0))+2);
977: scale = PetscPowRealInt(2.0,-sexp);
978: PetscCallCUBLAS(cublasXscal(cublasv2handle,ld2,&scale,d_As,one));
979: PetscCall(PetscLogGpuFlops(1.0*n*n));
980: } else sexp = 0;
982: /* Horner evaluation */
983: PetscCallCUBLAS(cublasXgemm(cublasv2handle,CUBLAS_OP_N,CUBLAS_OP_N,n,n,n,&sone,d_As,ld,d_As,ld,&szero,d_A2,ld));
984: PetscCall(PetscLogGpuFlops(2.0*n*n*n));
985: PetscCallCUDA(cudaMemset(d_Q,0,sizeof(PetscScalar)*ld2));
986: PetscCallCUDA(cudaMemset(d_P,0,sizeof(PetscScalar)*ld2));
987: PetscCall(set_diagonal(n,d_Q,ld,c[p]));
988: PetscCall(set_diagonal(n,d_P,ld,c[p-1]));
990: odd = PETSC_TRUE;
991: for (k=p-1;k>0;k--) {
992: if (odd) {
993: PetscCallCUBLAS(cublasXgemm(cublasv2handle,CUBLAS_OP_N,CUBLAS_OP_N,n,n,n,&sone,d_Q,ld,d_A2,ld,&szero,d_W,ld));
994: SlepcSwap(d_Q,d_W,aux);
995: PetscCall(shift_diagonal(n,d_Q,ld,c[k-1]));
996: odd = PETSC_FALSE;
997: } else {
998: PetscCallCUBLAS(cublasXgemm(cublasv2handle,CUBLAS_OP_N,CUBLAS_OP_N,n,n,n,&sone,d_P,ld,d_A2,ld,&szero,d_W,ld));
999: SlepcSwap(d_P,d_W,aux);
1000: PetscCall(shift_diagonal(n,d_P,ld,c[k-1]));
1001: odd = PETSC_TRUE;
1002: }
1003: PetscCall(PetscLogGpuFlops(2.0*n*n*n));
1004: }
1005: if (odd) {
1006: PetscCallCUBLAS(cublasXgemm(cublasv2handle,CUBLAS_OP_N,CUBLAS_OP_N,n,n,n,&sone,d_Q,ld,d_As,ld,&szero,d_W,ld));
1007: SlepcSwap(d_Q,d_W,aux);
1008: PetscCallCUBLAS(cublasXaxpy(cublasv2handle,ld2,&smone,d_P,one,d_Q,one));
1010: ppQ[0] = d_Q;
1011: ppP[0] = d_P;
1012: PetscCallCUDA(cudaMemcpy(d_ppQ,ppQ,sizeof(PetscScalar*),cudaMemcpyHostToDevice));
1013: PetscCallCUDA(cudaMemcpy(d_ppP,ppP,sizeof(PetscScalar*),cudaMemcpyHostToDevice));
1015: PetscCallCUBLAS(cublasXgetrfBatched(cublasv2handle,n,d_ppQ,ld,d_ipiv,d_info,one));
1016: PetscCallCUDA(cudaMemcpy(&info,d_info,sizeof(PetscBLASInt),cudaMemcpyDeviceToHost));
1017: PetscCheck(info>=0,PETSC_COMM_SELF,PETSC_ERR_LIB,"LAPACKgetrf: Illegal value on argument %" PetscBLASInt_FMT,PetscAbsInt(info));
1018: PetscCheck(info<=0,PETSC_COMM_SELF,PETSC_ERR_MAT_LU_ZRPVT,"LAPACKgetrf: Matrix is singular. U(%" PetscBLASInt_FMT ",%" PetscBLASInt_FMT ") is zero",info,info);
1019: PetscCallCUBLAS(cublasXgetrsBatched(cublasv2handle,CUBLAS_OP_N,n,n,(const PetscScalar **)d_ppQ,ld,d_ipiv,d_ppP,ld,&info,one));
1020: PetscCheck(info>=0,PETSC_COMM_SELF,PETSC_ERR_LIB,"LAPACKgetri: Illegal value on argument %" PetscBLASInt_FMT,PetscAbsInt(info));
1021: PetscCheck(info<=0,PETSC_COMM_SELF,PETSC_ERR_MAT_LU_ZRPVT,"LAPACKgetri: Matrix is singular. U(%" PetscBLASInt_FMT ",%" PetscBLASInt_FMT ") is zero",info,info);
1022: PetscCallCUBLAS(cublasXscal(cublasv2handle,ld2,&stwo,d_P,one));
1023: PetscCall(shift_diagonal(n,d_P,ld,sone));
1024: PetscCallCUBLAS(cublasXscal(cublasv2handle,ld2,&smone,d_P,one));
1025: } else {
1026: PetscCallCUBLAS(cublasXgemm(cublasv2handle,CUBLAS_OP_N,CUBLAS_OP_N,n,n,n,&sone,d_P,ld,d_As,ld,&szero,d_W,ld));
1027: SlepcSwap(d_P,d_W,aux);
1028: PetscCallCUBLAS(cublasXaxpy(cublasv2handle,ld2,&smone,d_P,one,d_Q,one));
1030: ppQ[0] = d_Q;
1031: ppP[0] = d_P;
1032: PetscCallCUDA(cudaMemcpy(d_ppQ,ppQ,sizeof(PetscScalar*),cudaMemcpyHostToDevice));
1033: PetscCallCUDA(cudaMemcpy(d_ppP,ppP,sizeof(PetscScalar*),cudaMemcpyHostToDevice));
1035: PetscCallCUBLAS(cublasXgetrfBatched(cublasv2handle,n,d_ppQ,ld,d_ipiv,d_info,one));
1036: PetscCallCUDA(cudaMemcpy(&info,d_info,sizeof(PetscBLASInt),cudaMemcpyDeviceToHost));
1037: PetscCheck(info>=0,PETSC_COMM_SELF,PETSC_ERR_LIB,"LAPACKgetrf: Illegal value on argument %" PetscBLASInt_FMT,PetscAbsInt(info));
1038: PetscCheck(info<=0,PETSC_COMM_SELF,PETSC_ERR_MAT_LU_ZRPVT,"LAPACKgetrf: Matrix is singular. U(%" PetscBLASInt_FMT ",%" PetscBLASInt_FMT ") is zero",info,info);
1039: PetscCallCUBLAS(cublasXgetrsBatched(cublasv2handle,CUBLAS_OP_N,n,n,(const PetscScalar **)d_ppQ,ld,d_ipiv,d_ppP,ld,&info,one));
1040: PetscCheck(info>=0,PETSC_COMM_SELF,PETSC_ERR_LIB,"LAPACKgetri: Illegal value on argument %" PetscBLASInt_FMT,PetscAbsInt(info));
1041: PetscCheck(info<=0,PETSC_COMM_SELF,PETSC_ERR_MAT_LU_ZRPVT,"LAPACKgetri: Matrix is singular. U(%" PetscBLASInt_FMT ",%" PetscBLASInt_FMT ") is zero",info,info);
1042: PetscCallCUBLAS(cublasXscal(cublasv2handle,ld2,&stwo,d_P,one));
1043: PetscCall(shift_diagonal(n,d_P,ld,sone));
1044: }
1045: PetscCall(PetscLogGpuFlops(2.0*n*n*n+2.0*n*n*n/3.0+4.0*n*n));
1047: for (k=1;k<=sexp;k++) {
1048: PetscCallCUBLAS(cublasXgemm(cublasv2handle,CUBLAS_OP_N,CUBLAS_OP_N,n,n,n,&sone,d_P,ld,d_P,ld,&szero,d_W,ld));
1049: PetscCallCUDA(cudaMemcpy(d_P,d_W,sizeof(PetscScalar)*ld2,cudaMemcpyDeviceToDevice));
1050: }
1051: PetscCall(PetscLogGpuFlops(2.0*n*n*n*sexp));
1053: PetscCall(PetscLogGpuTimeEnd());
1054: PetscCallCUDA(cudaFree(d_Q));
1055: PetscCallCUDA(cudaFree(d_W));
1056: PetscCallCUDA(cudaFree(d_A2));
1057: PetscCallCUDA(cudaFree(d_ipiv));
1058: PetscCallCUDA(cudaFree(d_info));
1059: PetscCallCUDA(cudaFree(d_ppP));
1060: PetscCallCUDA(cudaFree(d_ppQ));
1062: PetscCall(PetscFree(ppP));
1063: PetscCall(PetscFree(ppQ));
1065: if (d_P!=d_Ba) PetscCallCUDA(cudaMemcpy(d_Ba,d_P,sizeof(PetscScalar)*ld2,cudaMemcpyDeviceToDevice));
1066: if (A!=B) {
1067: if (s>0.5) { /* undo scaling */
1068: scale = 1.0/scale;
1069: PetscCallCUBLAS(cublasXscal(cublasv2handle,ld2,&scale,d_As,one));
1070: }
1071: PetscCall(MatDenseCUDARestoreArrayRead(A,(const PetscScalar**)&d_As));
1072: } else PetscCallCUDA(cudaFree(d_As));
1073: PetscCall(MatDenseCUDARestoreArrayWrite(B,&d_Ba));
1074: PetscFunctionReturn(PETSC_SUCCESS);
1075: }
1077: #if PetscDefined(HAVE_MAGMA)
1078: #include <slepcmagma.h>
1080: PetscErrorCode FNEvaluateFunctionMat_Exp_Pade_CUDAm(FN fn,Mat A,Mat B)
1081: {
1082: PetscBLASInt n=0,ld,ld2,*piv,one=1;
1083: PetscInt m,k,sexp;
1084: PetscBool odd;
1085: const PetscInt p=MAX_PADE;
1086: PetscReal c[MAX_PADE+1],s;
1087: PetscScalar scale,smone=-1.0,sone=1.0,stwo=2.0,szero=0.0;
1088: const PetscScalar *Aa;
1089: PetscScalar *d_Ba,*d_As,*d_A2,*d_Q,*d_P,*d_W,*aux;
1090: cublasHandle_t cublasv2handle;
1092: PetscFunctionBegin;
1093: PetscCall(PetscDeviceInitialize(PETSC_DEVICE_CUDA)); /* For CUDA event timers */
1094: PetscCall(PetscCUBLASGetHandle(&cublasv2handle));
1095: PetscCall(SlepcMagmaInit());
1096: PetscCall(MatGetSize(A,&m,NULL));
1097: PetscCall(PetscBLASIntCast(m,&n));
1098: ld = n;
1099: ld2 = ld*ld;
1100: if (A==B) {
1101: PetscCallCUDA(cudaMalloc((void **)&d_As,sizeof(PetscScalar)*m*m));
1102: PetscCall(MatDenseCUDAGetArrayRead(A,&Aa));
1103: PetscCallCUDA(cudaMemcpy(d_As,Aa,sizeof(PetscScalar)*ld2,cudaMemcpyDeviceToDevice));
1104: PetscCall(MatDenseCUDARestoreArrayRead(A,&Aa));
1105: } else PetscCall(MatDenseCUDAGetArrayRead(A,(const PetscScalar**)&d_As));
1106: PetscCall(MatDenseCUDAGetArrayWrite(B,&d_Ba));
1108: PetscCallCUDA(cudaMalloc((void **)&d_Q,sizeof(PetscScalar)*m*m));
1109: PetscCallCUDA(cudaMalloc((void **)&d_W,sizeof(PetscScalar)*m*m));
1110: PetscCallCUDA(cudaMalloc((void **)&d_A2,sizeof(PetscScalar)*m*m));
1112: PetscCall(PetscMalloc1(n,&piv));
1114: d_P = d_Ba;
1115: PetscCall(PetscLogGpuTimeBegin());
1117: /* Pade' coefficients */
1118: c[0] = 1.0;
1119: for (k=1;k<=p;k++) c[k] = c[k-1]*(p+1-k)/(k*(2*p+1-k));
1121: /* Scaling */
1122: PetscCallCUBLAS(cublasXnrm2(cublasv2handle,ld2,d_As,one,&s));
1123: PetscCall(PetscLogGpuFlops(1.0*n*n));
1125: if (s>0.5) {
1126: sexp = PetscMax(0,(int)(PetscLogReal(s)/PetscLogReal(2.0))+2);
1127: scale = PetscPowRealInt(2.0,-sexp);
1128: PetscCallCUBLAS(cublasXscal(cublasv2handle,ld2,&scale,d_As,one));
1129: PetscCall(PetscLogGpuFlops(1.0*n*n));
1130: } else sexp = 0;
1132: /* Horner evaluation */
1133: PetscCallCUBLAS(cublasXgemm(cublasv2handle,CUBLAS_OP_N,CUBLAS_OP_N,n,n,n,&sone,d_As,ld,d_As,ld,&szero,d_A2,ld));
1134: PetscCall(PetscLogGpuFlops(2.0*n*n*n));
1135: PetscCallCUDA(cudaMemset(d_Q,0,sizeof(PetscScalar)*ld2));
1136: PetscCallCUDA(cudaMemset(d_P,0,sizeof(PetscScalar)*ld2));
1137: PetscCall(set_diagonal(n,d_Q,ld,c[p]));
1138: PetscCall(set_diagonal(n,d_P,ld,c[p-1]));
1140: odd = PETSC_TRUE;
1141: for (k=p-1;k>0;k--) {
1142: if (odd) {
1143: PetscCallCUBLAS(cublasXgemm(cublasv2handle,CUBLAS_OP_N,CUBLAS_OP_N,n,n,n,&sone,d_Q,ld,d_A2,ld,&szero,d_W,ld));
1144: SlepcSwap(d_Q,d_W,aux);
1145: PetscCall(shift_diagonal(n,d_Q,ld,c[k-1]));
1146: odd = PETSC_FALSE;
1147: } else {
1148: PetscCallCUBLAS(cublasXgemm(cublasv2handle,CUBLAS_OP_N,CUBLAS_OP_N,n,n,n,&sone,d_P,ld,d_A2,ld,&szero,d_W,ld));
1149: SlepcSwap(d_P,d_W,aux);
1150: PetscCall(shift_diagonal(n,d_P,ld,c[k-1]));
1151: odd = PETSC_TRUE;
1152: }
1153: PetscCall(PetscLogGpuFlops(2.0*n*n*n));
1154: }
1155: if (odd) {
1156: PetscCallCUBLAS(cublasXgemm(cublasv2handle,CUBLAS_OP_N,CUBLAS_OP_N,n,n,n,&sone,d_Q,ld,d_As,ld,&szero,d_W,ld));
1157: SlepcSwap(d_Q,d_W,aux);
1158: PetscCallCUBLAS(cublasXaxpy(cublasv2handle,ld2,&smone,d_P,one,d_Q,one));
1159: PetscCallMAGMA(magma_xgesv_gpu,n,n,d_Q,ld,piv,d_P,ld);
1160: PetscCallCUBLAS(cublasXscal(cublasv2handle,ld2,&stwo,d_P,one));
1161: PetscCall(shift_diagonal(n,d_P,ld,sone));
1162: PetscCallCUBLAS(cublasXscal(cublasv2handle,ld2,&smone,d_P,one));
1163: } else {
1164: PetscCallCUBLAS(cublasXgemm(cublasv2handle,CUBLAS_OP_N,CUBLAS_OP_N,n,n,n,&sone,d_P,ld,d_As,ld,&szero,d_W,ld));
1165: SlepcSwap(d_P,d_W,aux);
1166: PetscCallCUBLAS(cublasXaxpy(cublasv2handle,ld2,&smone,d_P,one,d_Q,one));
1167: PetscCallMAGMA(magma_xgesv_gpu,n,n,d_Q,ld,piv,d_P,ld);
1168: PetscCallCUBLAS(cublasXscal(cublasv2handle,ld2,&stwo,d_P,one));
1169: PetscCall(shift_diagonal(n,d_P,ld,sone));
1170: }
1171: PetscCall(PetscLogGpuFlops(2.0*n*n*n+2.0*n*n*n/3.0+4.0*n*n));
1173: for (k=1;k<=sexp;k++) {
1174: PetscCallCUBLAS(cublasXgemm(cublasv2handle,CUBLAS_OP_N,CUBLAS_OP_N,n,n,n,&sone,d_P,ld,d_P,ld,&szero,d_W,ld));
1175: PetscCallCUDA(cudaMemcpy(d_P,d_W,sizeof(PetscScalar)*ld2,cudaMemcpyDeviceToDevice));
1176: }
1177: PetscCall(PetscLogGpuFlops(2.0*n*n*n*sexp));
1179: PetscCall(PetscLogGpuTimeEnd());
1180: PetscCallCUDA(cudaFree(d_Q));
1181: PetscCallCUDA(cudaFree(d_W));
1182: PetscCallCUDA(cudaFree(d_A2));
1183: PetscCall(PetscFree(piv));
1185: if (d_P!=d_Ba) PetscCallCUDA(cudaMemcpy(d_Ba,d_P,sizeof(PetscScalar)*ld2,cudaMemcpyDeviceToDevice));
1186: if (A!=B) {
1187: if (s>0.5) { /* undo scaling */
1188: scale = 1.0/scale;
1189: PetscCallCUBLAS(cublasXscal(cublasv2handle,ld2,&scale,d_As,one));
1190: }
1191: PetscCall(MatDenseCUDARestoreArrayRead(A,(const PetscScalar**)&d_As));
1192: } else PetscCallCUDA(cudaFree(d_As));
1193: PetscCall(MatDenseCUDARestoreArrayWrite(B,&d_Ba));
1194: PetscFunctionReturn(PETSC_SUCCESS);
1195: }
1197: /*
1198: * Matrix exponential implementation based on algorithm and matlab code by N. Higham and co-authors
1199: *
1200: * N. J. Higham, "The scaling and squaring method for the matrix exponential
1201: * revisited", SIAM J. Matrix Anal. Appl. 26(4):1179-1193, 2005.
1202: */
1203: PetscErrorCode FNEvaluateFunctionMat_Exp_Higham_CUDAm(FN fn,Mat A,Mat B)
1204: {
1205: PetscBLASInt n_=0,n2,*ipiv,one=1;
1206: PetscInt n,m,j,s;
1207: PetscScalar scale,smone=-1.0,sone=1.0,stwo=2.0,szero=0.0;
1208: PetscScalar *d_Ba,*Apowers[5],*d_Apowers[5],*d_Q,*d_P,*d_W,*work,*d_work,*aux;
1209: const PetscScalar *Aa,*c;
1210: const PetscScalar c3[4] = { 120, 60, 12, 1 };
1211: const PetscScalar c5[6] = { 30240, 15120, 3360, 420, 30, 1 };
1212: const PetscScalar c7[8] = { 17297280, 8648640, 1995840, 277200, 25200, 1512, 56, 1 };
1213: const PetscScalar c9[10] = { 17643225600, 8821612800, 2075673600, 302702400, 30270240,
1214: 2162160, 110880, 3960, 90, 1 };
1215: const PetscScalar c13[14] = { 64764752532480000, 32382376266240000, 7771770303897600,
1216: 1187353796428800, 129060195264000, 10559470521600,
1217: 670442572800, 33522128640, 1323241920,
1218: 40840800, 960960, 16380, 182, 1 };
1219: cublasHandle_t cublasv2handle;
1221: PetscFunctionBegin;
1222: PetscCall(PetscDeviceInitialize(PETSC_DEVICE_CUDA)); /* For CUDA event timers */
1223: PetscCall(PetscCUBLASGetHandle(&cublasv2handle));
1224: PetscCall(SlepcMagmaInit());
1225: PetscCall(MatGetSize(A,&n,NULL));
1226: PetscCall(PetscBLASIntCast(n,&n_));
1227: n2 = n_*n_;
1228: PetscCall(PetscMalloc2(8*n*n,&work,n,&ipiv));
1229: /* Matrix powers */
1230: Apowers[0] = work; /* Apowers[0] = A */
1231: Apowers[1] = Apowers[0] + n*n; /* Apowers[1] = A^2 */
1232: Apowers[2] = Apowers[1] + n*n; /* Apowers[2] = A^4 */
1233: Apowers[3] = Apowers[2] + n*n; /* Apowers[3] = A^6 */
1234: Apowers[4] = Apowers[3] + n*n; /* Apowers[4] = A^8 */
1235: if (A==B) {
1236: PetscCallCUDA(cudaMalloc((void**)&d_work,7*n*n*sizeof(PetscScalar)));
1237: d_Apowers[0] = d_work; /* d_Apowers[0] = A */
1238: d_Apowers[1] = d_Apowers[0] + n*n; /* d_Apowers[1] = A^2 */
1239: PetscCall(MatDenseCUDAGetArrayRead(A,&Aa));
1240: PetscCallCUDA(cudaMemcpy(d_Apowers[0],Aa,n2*sizeof(PetscScalar),cudaMemcpyDeviceToDevice));
1241: PetscCall(MatDenseCUDARestoreArrayRead(A,&Aa));
1242: } else {
1243: PetscCallCUDA(cudaMalloc((void**)&d_work,6*n*n*sizeof(PetscScalar)));
1244: PetscCall(MatDenseCUDAGetArrayRead(A,(const PetscScalar**)&d_Apowers[0]));
1245: d_Apowers[1] = d_work; /* d_Apowers[1] = A^2 */
1246: }
1247: PetscCall(MatDenseCUDAGetArrayWrite(B,&d_Ba));
1248: d_Apowers[2] = d_Apowers[1] + n*n; /* d_Apowers[2] = A^4 */
1249: d_Apowers[3] = d_Apowers[2] + n*n; /* d_Apowers[3] = A^6 */
1250: d_Apowers[4] = d_Apowers[3] + n*n; /* d_Apowers[4] = A^8 */
1251: d_Q = d_Apowers[4] + n*n;
1252: d_W = d_Q + n*n;
1254: PetscCall(PetscLogGpuTimeBegin());
1256: PetscCallCUBLAS(cublasXgemm(cublasv2handle,CUBLAS_OP_N,CUBLAS_OP_N,n_,n_,n_,&sone,d_Apowers[0],n_,d_Apowers[0],n_,&szero,d_Apowers[1],n_));
1257: PetscCallCUBLAS(cublasXgemm(cublasv2handle,CUBLAS_OP_N,CUBLAS_OP_N,n_,n_,n_,&sone,d_Apowers[1],n_,d_Apowers[1],n_,&szero,d_Apowers[2],n_));
1258: PetscCallCUBLAS(cublasXgemm(cublasv2handle,CUBLAS_OP_N,CUBLAS_OP_N,n_,n_,n_,&sone,d_Apowers[1],n_,d_Apowers[2],n_,&szero,d_Apowers[3],n_));
1259: PetscCall(PetscLogGpuFlops(6.0*n*n*n));
1261: PetscCallCUDA(cudaMemcpy(Apowers[0],d_Apowers[0],n2*sizeof(PetscScalar),cudaMemcpyDeviceToHost));
1262: PetscCallCUDA(cudaMemcpy(Apowers[1],d_Apowers[1],3*n2*sizeof(PetscScalar),cudaMemcpyDeviceToHost));
1263: PetscCall(PetscLogGpuToCpu(4*n2*sizeof(PetscScalar)));
1264: /* Compute scaling parameter and order of Pade approximant */
1265: PetscCall(expm_params(n,Apowers,&s,&m,Apowers[4]));
1267: if (s) { /* rescale */
1268: for (j=0;j<4;j++) {
1269: scale = PetscPowRealInt(2.0,-PetscMax(2*j,1)*s);
1270: PetscCallCUBLAS(cublasXscal(cublasv2handle,n2,&scale,d_Apowers[j],one));
1271: }
1272: PetscCall(PetscLogGpuFlops(4.0*n*n));
1273: }
1275: /* Evaluate the Pade approximant */
1276: switch (m) {
1277: case 3: c = c3; break;
1278: case 5: c = c5; break;
1279: case 7: c = c7; break;
1280: case 9: c = c9; break;
1281: case 13: c = c13; break;
1282: default: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Wrong value of m %" PetscInt_FMT,m);
1283: }
1284: d_P = d_Ba;
1285: switch (m) {
1286: case 3:
1287: case 5:
1288: case 7:
1289: case 9:
1290: if (m==9) PetscCallCUBLAS(cublasXgemm(cublasv2handle,CUBLAS_OP_N,CUBLAS_OP_N,n_,n_,n_,&sone,d_Apowers[1],n_,d_Apowers[3],n_,&szero,d_Apowers[4],n_));
1291: PetscCallCUDA(cudaMemset(d_P,0,sizeof(PetscScalar)*n2));
1292: PetscCallCUDA(cudaMemset(d_Q,0,sizeof(PetscScalar)*n2));
1293: PetscCall(set_diagonal(n,d_P,n,c[1]));
1294: PetscCall(set_diagonal(n,d_Q,n,c[0]));
1295: for (j=m;j>=3;j-=2) {
1296: PetscCallCUBLAS(cublasXaxpy(cublasv2handle,n2,&c[j],d_Apowers[(j+1)/2-1],one,d_P,one));
1297: PetscCallCUBLAS(cublasXaxpy(cublasv2handle,n2,&c[j-1],d_Apowers[(j+1)/2-1],one,d_Q,one));
1298: PetscCall(PetscLogGpuFlops(4.0*n*n));
1299: }
1300: PetscCallCUBLAS(cublasXgemm(cublasv2handle,CUBLAS_OP_N,CUBLAS_OP_N,n_,n_,n_,&sone,d_Apowers[0],n_,d_P,n_,&szero,d_W,n_));
1301: PetscCall(PetscLogGpuFlops(2.0*n*n*n));
1302: SlepcSwap(d_P,d_W,aux);
1303: break;
1304: case 13:
1305: /* P = A*(Apowers[3]*(c[13]*Apowers[3] + c[11]*Apowers[2] + c[9]*Apowers[1])
1306: + c[7]*Apowers[3] + c[5]*Apowers[2] + c[3]*Apowers[1] + c[1]*I) */
1307: PetscCallCUDA(cudaMemcpy(d_P,d_Apowers[3],n2*sizeof(PetscScalar),cudaMemcpyDeviceToDevice));
1308: PetscCallCUBLAS(cublasXscal(cublasv2handle,n2,&c[13],d_P,one));
1309: PetscCallCUBLAS(cublasXaxpy(cublasv2handle,n2,&c[11],d_Apowers[2],one,d_P,one));
1310: PetscCallCUBLAS(cublasXaxpy(cublasv2handle,n2,&c[9],d_Apowers[1],one,d_P,one));
1311: PetscCallCUBLAS(cublasXgemm(cublasv2handle,CUBLAS_OP_N,CUBLAS_OP_N,n_,n_,n_,&sone,d_Apowers[3],n_,d_P,n_,&szero,d_W,n_));
1312: PetscCall(PetscLogGpuFlops(5.0*n*n+2.0*n*n*n));
1314: PetscCallCUDA(cudaMemset(d_P,0,sizeof(PetscScalar)*n2));
1315: PetscCall(set_diagonal(n,d_P,n,c[1]));
1316: PetscCallCUBLAS(cublasXaxpy(cublasv2handle,n2,&c[7],d_Apowers[3],one,d_P,one));
1317: PetscCallCUBLAS(cublasXaxpy(cublasv2handle,n2,&c[5],d_Apowers[2],one,d_P,one));
1318: PetscCallCUBLAS(cublasXaxpy(cublasv2handle,n2,&c[3],d_Apowers[1],one,d_P,one));
1319: PetscCallCUBLAS(cublasXaxpy(cublasv2handle,n2,&sone,d_P,one,d_W,one));
1320: PetscCallCUBLAS(cublasXgemm(cublasv2handle,CUBLAS_OP_N,CUBLAS_OP_N,n_,n_,n_,&sone,d_Apowers[0],n_,d_W,n_,&szero,d_P,n_));
1321: PetscCall(PetscLogGpuFlops(7.0*n*n+2.0*n*n*n));
1322: /* Q = Apowers[3]*(c[12]*Apowers[3] + c[10]*Apowers[2] + c[8]*Apowers[1])
1323: + c[6]*Apowers[3] + c[4]*Apowers[2] + c[2]*Apowers[1] + c[0]*I */
1324: PetscCallCUDA(cudaMemcpy(d_Q,d_Apowers[3],n2*sizeof(PetscScalar),cudaMemcpyDeviceToDevice));
1325: PetscCallCUBLAS(cublasXscal(cublasv2handle,n2,&c[12],d_Q,one));
1326: PetscCallCUBLAS(cublasXaxpy(cublasv2handle,n2,&c[10],d_Apowers[2],one,d_Q,one));
1327: PetscCallCUBLAS(cublasXaxpy(cublasv2handle,n2,&c[8],d_Apowers[1],one,d_Q,one));
1328: PetscCallCUBLAS(cublasXgemm(cublasv2handle,CUBLAS_OP_N,CUBLAS_OP_N,n_,n_,n_,&sone,d_Apowers[3],n_,d_Q,n_,&szero,d_W,n_));
1329: PetscCall(PetscLogGpuFlops(5.0*n*n+2.0*n*n*n));
1330: PetscCallCUDA(cudaMemset(d_Q,0,sizeof(PetscScalar)*n2));
1331: PetscCall(set_diagonal(n,d_Q,n,c[0]));
1332: PetscCallCUBLAS(cublasXaxpy(cublasv2handle,n2,&c[6],d_Apowers[3],one,d_Q,one));
1333: PetscCallCUBLAS(cublasXaxpy(cublasv2handle,n2,&c[4],d_Apowers[2],one,d_Q,one));
1334: PetscCallCUBLAS(cublasXaxpy(cublasv2handle,n2,&c[2],d_Apowers[1],one,d_Q,one));
1335: PetscCallCUBLAS(cublasXaxpy(cublasv2handle,n2,&sone,d_W,one,d_Q,one));
1336: PetscCall(PetscLogGpuFlops(7.0*n*n));
1337: break;
1338: default: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Wrong value of m %" PetscInt_FMT,m);
1339: }
1340: PetscCallCUBLAS(cublasXaxpy(cublasv2handle,n2,&smone,d_P,one,d_Q,one));
1342: PetscCallMAGMA(magma_xgesv_gpu,n_,n_,d_Q,n_,ipiv,d_P,n_);
1344: PetscCallCUBLAS(cublasXscal(cublasv2handle,n2,&stwo,d_P,one));
1345: PetscCall(shift_diagonal(n,d_P,n,sone));
1346: PetscCall(PetscLogGpuFlops(2.0*n*n*n/3.0+4.0*n*n));
1348: /* Squaring */
1349: for (j=1;j<=s;j++) {
1350: PetscCallCUBLAS(cublasXgemm(cublasv2handle,CUBLAS_OP_N,CUBLAS_OP_N,n_,n_,n_,&sone,d_P,n_,d_P,n_,&szero,d_W,n_));
1351: SlepcSwap(d_P,d_W,aux);
1352: }
1353: PetscCall(PetscLogGpuFlops(2.0*n*n*n*s));
1354: PetscCall(PetscLogGpuTimeEnd());
1356: PetscCall(PetscFree2(work,ipiv));
1357: if (d_P!=d_Ba) PetscCallCUDA(cudaMemcpy(d_Ba,d_P,n2*sizeof(PetscScalar),cudaMemcpyDeviceToDevice));
1358: if (A!=B) {
1359: if (s>0.5) { /* undo scaling */
1360: scale = 1.0/PetscPowRealInt(2.0,-s);
1361: PetscCallCUBLAS(cublasXscal(cublasv2handle,n2,&scale,d_Apowers[0],one));
1362: }
1363: PetscCall(MatDenseCUDARestoreArrayRead(A,(const PetscScalar**)&d_Apowers[0]));
1364: }
1365: PetscCall(MatDenseCUDARestoreArrayWrite(B,&d_Ba));
1366: PetscCallCUDA(cudaFree(d_work));
1367: PetscFunctionReturn(PETSC_SUCCESS);
1368: }
1370: /*
1371: * Matrix exponential implementation based on algorithm and matlab code by Stefan Guettel
1372: * and Yuji Nakatsukasa
1373: *
1374: * Stefan Guettel and Yuji Nakatsukasa, "Scaled and Squared Subdiagonal Pade'
1375: * Approximation for the Matrix Exponential",
1376: * SIAM J. Matrix Anal. Appl. 37(1):145-170, 2016.
1377: * https://doi.org/10.1137/15M1027553
1378: */
1379: PetscErrorCode FNEvaluateFunctionMat_Exp_GuettelNakatsukasa_CUDAm(FN fn,Mat A,Mat B)
1380: {
1381: PetscInt i,j,n_,s,k,m,mod;
1382: PetscBLASInt n=0,n2=0,irsize=0,rsizediv2,ipsize=0,iremainsize=0,query=-1,*piv,minlen,lwork=0,one=1;
1383: PetscReal nrm,shift=0.0,rone=1.0,rzero=0.0;
1384: #if PetscDefined(USE_COMPLEX)
1385: PetscReal *rwork=NULL;
1386: #endif
1387: PetscComplex *d_As,*d_RR,*d_RR2,*d_expmA,*d_expmA2,*d_Maux,*d_Maux2,rsize,*r,psize,*p,remainsize,*remainterm,*rootp,*rootq,mult=0.0,scale,cone=1.0,czero=0.0,*aux;
1388: PetscScalar *d_Aa,*d_Ba,*d_Ba2,*Maux,*d_sMaux,*wr,*wi,expshift,sone=1.0,szero=0.0,*work,work1,*saux;
1389: const PetscScalar *Aa;
1390: PetscBool isreal,*d_isreal,flg;
1391: cublasHandle_t cublasv2handle;
1393: PetscFunctionBegin;
1394: PetscCall(PetscDeviceInitialize(PETSC_DEVICE_CUDA)); /* For CUDA event timers */
1395: PetscCall(PetscCUBLASGetHandle(&cublasv2handle));
1396: PetscCall(SlepcMagmaInit());
1397: PetscCall(MatGetSize(A,&n_,NULL));
1398: PetscCall(PetscBLASIntCast(n_,&n));
1399: PetscCall(PetscBLASIntCast(n*n,&n2));
1401: if (A==B) {
1402: PetscCallCUDA(cudaMalloc((void **)&d_Aa,sizeof(PetscScalar)*n2));
1403: PetscCall(MatDenseCUDAGetArrayRead(A,&Aa));
1404: PetscCallCUDA(cudaMemcpy(d_Aa,Aa,sizeof(PetscScalar)*n2,cudaMemcpyDeviceToDevice));
1405: PetscCall(MatDenseCUDARestoreArrayRead(A,&Aa));
1406: } else PetscCall(MatDenseCUDAGetArrayRead(A,(const PetscScalar**)&d_Aa));
1407: PetscCall(MatDenseCUDAGetArrayWrite(B,&d_Ba));
1408: d_Ba2 = d_Ba;
1410: PetscCallCUDA(cudaMalloc((void **)&d_isreal,sizeof(PetscBool)));
1411: PetscCallCUDA(cudaMalloc((void **)&d_sMaux,sizeof(PetscScalar)*n2));
1412: PetscCallCUDA(cudaMalloc((void **)&d_Maux,sizeof(PetscComplex)*n2));
1414: PetscCall(PetscLogGpuTimeBegin());
1415: d_Maux2 = d_Maux;
1416: PetscCall(PetscOptionsGetReal(NULL,NULL,"-fn_expm_estimated_eig",&shift,&flg));
1417: if (!flg) {
1418: PetscCall(PetscMalloc2(n,&wr,n,&wi));
1419: /* estimate rightmost eigenvalue and shift A with it */
1420: PetscCall(PetscMalloc1(n2,&Maux));
1421: PetscCall(MatDenseGetArrayRead(A,&Aa));
1422: PetscCall(PetscArraycpy(Maux,Aa,n2));
1423: PetscCall(MatDenseRestoreArrayRead(A,&Aa));
1424: #if !PetscDefined(USE_COMPLEX)
1425: PetscCallMAGMA(magma_xgeev,MagmaNoVec,MagmaNoVec,n,Maux,n,wr,wi,NULL,n,NULL,n,&work1,query);
1426: PetscCall(PetscBLASIntCast((PetscInt)PetscRealPart(work1),&lwork));
1427: PetscCall(PetscMalloc1(lwork,&work));
1428: PetscCallMAGMA(magma_xgeev,MagmaNoVec,MagmaNoVec,n,Maux,n,wr,wi,NULL,n,NULL,n,work,lwork);
1429: PetscCall(PetscFree(work));
1430: #else
1431: PetscCallMAGMA(magma_xgeev,MagmaNoVec,MagmaNoVec,n,Maux,n,wr,NULL,n,NULL,n,&work1,query,rwork);
1432: PetscCall(PetscBLASIntCast((PetscInt)PetscRealPart(work1),&lwork));
1433: PetscCall(PetscMalloc2(2*n,&rwork,lwork,&work));
1434: PetscCallMAGMA(magma_xgeev,MagmaNoVec,MagmaNoVec,n,Maux,n,wr,NULL,n,NULL,n,work,lwork,rwork);
1435: PetscCall(PetscFree2(rwork,work));
1436: #endif
1437: PetscCall(PetscFree(Maux));
1438: PetscCall(PetscLogGpuFlops(25.0*n*n*n+(n*n*n)/3.0+1.0*n*n*n));
1440: shift = PetscRealPart(wr[0]);
1441: for (i=1;i<n;i++) {
1442: if (PetscRealPart(wr[i]) > shift) shift = PetscRealPart(wr[i]);
1443: }
1444: PetscCall(PetscFree2(wr,wi));
1445: }
1446: /* shift so that largest real part is (about) 0 */
1447: PetscCallCUDA(cudaMemcpy(d_sMaux,d_Aa,sizeof(PetscScalar)*n2,cudaMemcpyDeviceToDevice));
1448: if (shift) {
1449: PetscCall(shift_diagonal(n,d_sMaux,n,-shift));
1450: PetscCall(PetscLogGpuFlops(1.0*n));
1451: }
1452: #if PetscDefined(USE_COMPLEX)
1453: PetscCallCUDA(cudaMemcpy(d_Maux,d_Aa,sizeof(PetscScalar)*n2,cudaMemcpyDeviceToDevice));
1454: if (shift) {
1455: PetscCall(shift_diagonal(n,d_Maux,n,-shift));
1456: PetscCall(PetscLogGpuFlops(1.0*n));
1457: }
1458: #endif
1459: if (A!=B) PetscCall(MatDenseCUDARestoreArrayRead(A,(const PetscScalar**)&d_Aa));
1460: else PetscCallCUDA(cudaFree(d_Aa));
1462: /* estimate norm(A) and select the scaling factor */
1463: PetscCallCUBLAS(cublasXnrm2(cublasv2handle,n2,d_sMaux,one,&nrm));
1464: PetscCall(PetscLogGpuFlops(2.0*n*n));
1465: PetscCall(sexpm_params(nrm,&s,&k,&m));
1466: if (s==0 && k==1 && m==0) { /* exp(A) = I+A to eps! */
1467: if (shift) expshift = PetscExpReal(shift);
1468: PetscCall(shift_Cdiagonal(n,d_Maux,n,rone,rzero));
1469: if (shift) {
1470: PetscCallCUBLAS(cublasXscal(cublasv2handle,n2,&expshift,d_sMaux,one));
1471: PetscCall(PetscLogGpuFlops(1.0*(n+n2)));
1472: } else PetscCall(PetscLogGpuFlops(1.0*n));
1473: PetscCallCUDA(cudaMemcpy(d_Ba,d_sMaux,sizeof(PetscScalar)*n2,cudaMemcpyDeviceToDevice));
1474: PetscCallCUDA(cudaFree(d_isreal));
1475: PetscCallCUDA(cudaFree(d_sMaux));
1476: PetscCallCUDA(cudaFree(d_Maux));
1477: PetscCall(MatDenseCUDARestoreArrayWrite(B,&d_Ba));
1478: PetscFunctionReturn(PETSC_SUCCESS); /* quick return */
1479: }
1481: PetscCallCUDA(cudaMalloc((void **)&d_expmA,sizeof(PetscComplex)*n2));
1482: PetscCallCUDA(cudaMalloc((void **)&d_As,sizeof(PetscComplex)*n2));
1483: PetscCallCUDA(cudaMalloc((void **)&d_RR,sizeof(PetscComplex)*n2));
1484: d_expmA2 = d_expmA; d_RR2 = d_RR;
1485: PetscCall(PetscMalloc1(n,&piv));
1486: /* scale matrix */
1487: #if !PetscDefined(USE_COMPLEX)
1488: PetscCall(copy_array2D_S2C(n,n,d_As,n,d_sMaux,n));
1489: #else
1490: PetscCallCUDA(cudaMemcpy(d_As,d_sMaux,sizeof(PetscScalar)*n2,cudaMemcpyDeviceToDevice));
1491: #endif
1492: scale = 1.0/PetscPowRealInt(2.0,s);
1493: PetscCallCUBLAS(cublasXCscal(cublasv2handle,n2,(const cuComplex *)&scale,(cuComplex *)d_As,one));
1494: PetscCall(SlepcLogGpuFlopsComplex(1.0*n2));
1496: /* evaluate Pade approximant (partial fraction or product form) */
1497: if (fn->method==3 || !m) { /* partial fraction */
1498: PetscCall(getcoeffs(k,m,&rsize,&psize,&remainsize,PETSC_TRUE));
1499: PetscCall(PetscBLASIntCast((PetscInt)PetscRealPartComplex(rsize),&irsize));
1500: PetscCall(PetscBLASIntCast((PetscInt)PetscRealPartComplex(psize),&ipsize));
1501: PetscCall(PetscBLASIntCast((PetscInt)PetscRealPartComplex(remainsize),&iremainsize));
1502: PetscCall(PetscMalloc3(irsize,&r,ipsize,&p,iremainsize,&remainterm));
1503: PetscCall(getcoeffs(k,m,r,p,remainterm,PETSC_FALSE));
1505: PetscCallCUDA(cudaMemset(d_expmA,0,sizeof(PetscComplex)*n2));
1506: #if !PetscDefined(USE_COMPLEX)
1507: isreal = PETSC_TRUE;
1508: #else
1509: PetscCall(getisreal_array2D(n,n,d_Maux,n,d_isreal));
1510: PetscCallCUDA(cudaMemcpy(&isreal,d_isreal,sizeof(PetscBool),cudaMemcpyDeviceToHost));
1511: #endif
1512: if (isreal) {
1513: rsizediv2 = irsize/2;
1514: for (i=0;i<rsizediv2;i++) { /* use partial fraction to get R(As) */
1515: PetscCallCUDA(cudaMemcpy(d_Maux,d_As,sizeof(PetscComplex)*n2,cudaMemcpyDeviceToDevice));
1516: PetscCallCUDA(cudaMemset(d_RR,0,sizeof(PetscComplex)*n2));
1517: PetscCall(shift_Cdiagonal(n,d_Maux,n,-PetscRealPartComplex(p[2*i]),-PetscImaginaryPartComplex(p[2*i])));
1518: PetscCall(set_Cdiagonal(n,d_RR,n,PetscRealPartComplex(r[2*i]),PetscImaginaryPartComplex(r[2*i])));
1519: PetscCallMAGMA(magma_Cgesv_gpu,n,n,d_Maux,n,piv,d_RR,n);
1520: PetscCall(add_array2D_Conj(n,n,d_RR,n));
1521: PetscCallCUBLAS(cublasXCaxpy(cublasv2handle,n2,&cone,d_RR,one,d_expmA,one));
1522: /* shift(n) + gesv + axpy(n2) */
1523: PetscCall(SlepcLogGpuFlopsComplex(1.0*n+(2.0*n*n*n/3.0+2.0*n*n*n)+2.0*n2));
1524: }
1526: mod = ipsize % 2;
1527: if (mod) {
1528: PetscCallCUDA(cudaMemcpy(d_Maux,d_As,sizeof(PetscComplex)*n2,cudaMemcpyDeviceToDevice));
1529: PetscCallCUDA(cudaMemset(d_RR,0,sizeof(PetscComplex)*n2));
1530: PetscCall(shift_Cdiagonal(n,d_Maux,n,-PetscRealPartComplex(p[ipsize-1]),-PetscImaginaryPartComplex(p[ipsize-1])));
1531: PetscCall(set_Cdiagonal(n,d_RR,n,PetscRealPartComplex(r[irsize-1]),PetscImaginaryPartComplex(r[irsize-1])));
1532: PetscCallMAGMA(magma_Cgesv_gpu,n,n,d_Maux,n,piv,d_RR,n);
1533: PetscCallCUBLAS(cublasXCaxpy(cublasv2handle,n2,&cone,d_RR,one,d_expmA,one));
1534: PetscCall(SlepcLogGpuFlopsComplex(1.0*n+(2.0*n*n*n/3.0+2.0*n*n*n)+1.0*n2));
1535: }
1536: } else { /* complex */
1537: for (i=0;i<irsize;i++) { /* use partial fraction to get R(As) */
1538: PetscCallCUDA(cudaMemcpy(d_Maux,d_As,sizeof(PetscComplex)*n2,cudaMemcpyDeviceToDevice));
1539: PetscCallCUDA(cudaMemset(d_RR,0,sizeof(PetscComplex)*n2));
1540: PetscCall(shift_Cdiagonal(n,d_Maux,n,-PetscRealPartComplex(p[i]),-PetscImaginaryPartComplex(p[i])));
1541: PetscCall(set_Cdiagonal(n,d_RR,n,PetscRealPartComplex(r[i]),PetscImaginaryPartComplex(r[i])));
1542: PetscCallMAGMA(magma_Cgesv_gpu,n,n,d_Maux,n,piv,d_RR,n);
1543: PetscCallCUBLAS(cublasXCaxpy(cublasv2handle,n2,&cone,d_RR,one,d_expmA,one));
1544: PetscCall(SlepcLogGpuFlopsComplex(1.0*n+(2.0*n*n*n/3.0+2.0*n*n*n)+1.0*n2));
1545: }
1546: }
1547: for (i=0;i<iremainsize;i++) {
1548: if (!i) {
1549: PetscCallCUDA(cudaMemset(d_RR,0,sizeof(PetscComplex)*n2));
1550: PetscCall(set_Cdiagonal(n,d_RR,n,PetscRealPartComplex(remainterm[iremainsize-1]),PetscImaginaryPartComplex(remainterm[iremainsize-1])));
1551: } else {
1552: PetscCallCUDA(cudaMemcpy(d_RR,d_As,sizeof(PetscComplex)*n2,cudaMemcpyDeviceToDevice));
1553: for (j=1;j<i;j++) {
1554: PetscCallCUBLAS(cublasXCgemm(cublasv2handle,CUBLAS_OP_N,CUBLAS_OP_N,n,n,n,&cone,d_RR,n,d_RR,n,&czero,d_Maux,n));
1555: SlepcSwap(d_RR,d_Maux,aux);
1556: PetscCall(SlepcLogGpuFlopsComplex(2.0*n*n*n));
1557: }
1558: PetscCallCUBLAS(cublasXCscal(cublasv2handle,n2,&remainterm[iremainsize-1-i],d_RR,one));
1559: PetscCall(SlepcLogGpuFlopsComplex(1.0*n2));
1560: }
1561: PetscCallCUBLAS(cublasXCaxpy(cublasv2handle,n2,&cone,d_RR,one,d_expmA,one));
1562: PetscCall(SlepcLogGpuFlopsComplex(1.0*n2));
1563: }
1564: PetscCall(PetscFree3(r,p,remainterm));
1565: } else { /* product form, default */
1566: PetscCall(getcoeffsproduct(k,m,&rsize,&psize,&mult,PETSC_TRUE));
1567: PetscCall(PetscBLASIntCast((PetscInt)PetscRealPartComplex(rsize),&irsize));
1568: PetscCall(PetscBLASIntCast((PetscInt)PetscRealPartComplex(psize),&ipsize));
1569: PetscCall(PetscMalloc2(irsize,&rootp,ipsize,&rootq));
1570: PetscCall(getcoeffsproduct(k,m,rootp,rootq,&mult,PETSC_FALSE));
1572: PetscCallCUDA(cudaMemset(d_expmA,0,sizeof(PetscComplex)*n2));
1573: PetscCall(set_Cdiagonal(n,d_expmA,n,rone,rzero)); /* initialize */
1574: minlen = PetscMin(irsize,ipsize);
1575: for (i=0;i<minlen;i++) {
1576: PetscCallCUDA(cudaMemcpy(d_RR,d_As,sizeof(PetscComplex)*n2,cudaMemcpyDeviceToDevice));
1577: PetscCall(shift_Cdiagonal(n,d_RR,n,-PetscRealPartComplex(rootp[i]),-PetscImaginaryPartComplex(rootp[i])));
1578: PetscCallCUBLAS(cublasXCgemm(cublasv2handle,CUBLAS_OP_N,CUBLAS_OP_N,n,n,n,&cone,d_RR,n,d_expmA,n,&czero,d_Maux,n));
1579: SlepcSwap(d_expmA,d_Maux,aux);
1580: PetscCallCUDA(cudaMemcpy(d_RR,d_As,sizeof(PetscComplex)*n2,cudaMemcpyDeviceToDevice));
1581: PetscCall(shift_Cdiagonal(n,d_RR,n,-PetscRealPartComplex(rootq[i]),-PetscImaginaryPartComplex(rootq[i])));
1582: PetscCallMAGMA(magma_Cgesv_gpu,n,n,d_RR,n,piv,d_expmA,n);
1583: /* shift(n) + gemm + shift(n) + gesv */
1584: PetscCall(SlepcLogGpuFlopsComplex(1.0*n+(2.0*n*n*n)+1.0*n+(2.0*n*n*n/3.0+2.0*n*n*n)));
1585: }
1586: /* extra enumerator */
1587: for (i=minlen;i<irsize;i++) {
1588: PetscCallCUDA(cudaMemcpy(d_RR,d_As,sizeof(PetscComplex)*n2,cudaMemcpyDeviceToDevice));
1589: PetscCall(shift_Cdiagonal(n,d_RR,n,-PetscRealPartComplex(rootp[i]),-PetscImaginaryPartComplex(rootp[i])));
1590: PetscCallCUBLAS(cublasXCgemm(cublasv2handle,CUBLAS_OP_N,CUBLAS_OP_N,n,n,n,&cone,d_RR,n,d_expmA,n,&czero,d_Maux,n));
1591: SlepcSwap(d_expmA,d_Maux,aux);
1592: PetscCall(SlepcLogGpuFlopsComplex(1.0*n+2.0*n*n*n));
1593: }
1594: /* extra denominator */
1595: for (i=minlen;i<ipsize;i++) {
1596: PetscCallCUDA(cudaMemcpy(d_RR,d_As,sizeof(PetscComplex)*n2,cudaMemcpyDeviceToDevice));
1597: PetscCall(shift_Cdiagonal(n,d_RR,n,-PetscRealPartComplex(rootq[i]),-PetscImaginaryPartComplex(rootq[i])));
1598: PetscCallMAGMA(magma_Cgesv_gpu,n,n,d_RR,n,piv,d_expmA,n);
1599: PetscCall(SlepcLogGpuFlopsComplex(1.0*n+(2.0*n*n*n/3.0+2.0*n*n*n)));
1600: }
1601: PetscCallCUBLAS(cublasXCscal(cublasv2handle,n2,&mult,d_expmA,one));
1602: PetscCall(SlepcLogGpuFlopsComplex(1.0*n2));
1603: PetscCall(PetscFree2(rootp,rootq));
1604: }
1606: #if !PetscDefined(USE_COMPLEX)
1607: PetscCall(copy_array2D_C2S(n,n,d_Ba2,n,d_expmA,n));
1608: #else
1609: PetscCallCUDA(cudaMemcpy(d_Ba2,d_expmA,sizeof(PetscScalar)*n2,cudaMemcpyDeviceToDevice));
1610: #endif
1612: /* perform repeated squaring */
1613: for (i=0;i<s;i++) { /* final squaring */
1614: PetscCallCUBLAS(cublasXgemm(cublasv2handle,CUBLAS_OP_N,CUBLAS_OP_N,n,n,n,&sone,d_Ba2,n,d_Ba2,n,&szero,d_sMaux,n));
1615: SlepcSwap(d_Ba2,d_sMaux,saux);
1616: PetscCall(PetscLogGpuFlops(2.0*n*n*n));
1617: }
1618: if (d_Ba2!=d_Ba) {
1619: PetscCallCUDA(cudaMemcpy(d_Ba,d_Ba2,sizeof(PetscScalar)*n2,cudaMemcpyDeviceToDevice));
1620: d_sMaux = d_Ba2;
1621: }
1622: if (shift) {
1623: expshift = PetscExpReal(shift);
1624: PetscCallCUBLAS(cublasXscal(cublasv2handle,n2,&expshift,d_Ba,one));
1625: PetscCall(PetscLogGpuFlops(1.0*n2));
1626: }
1628: PetscCall(PetscLogGpuTimeEnd());
1630: /* restore pointers */
1631: d_Maux = d_Maux2; d_expmA = d_expmA2; d_RR = d_RR2;
1632: PetscCall(MatDenseCUDARestoreArrayWrite(B,&d_Ba));
1633: PetscCallCUDA(cudaFree(d_isreal));
1634: PetscCallCUDA(cudaFree(d_sMaux));
1635: PetscCallCUDA(cudaFree(d_Maux));
1636: PetscCallCUDA(cudaFree(d_expmA));
1637: PetscCallCUDA(cudaFree(d_As));
1638: PetscCallCUDA(cudaFree(d_RR));
1639: PetscCall(PetscFree(piv));
1640: PetscFunctionReturn(PETSC_SUCCESS);
1641: }
1642: #endif /* PETSC_HAVE_MAGMA */
1643: #endif /* PETSC_HAVE_CUDA */
1645: static PetscErrorCode FNView_Exp(FN fn,PetscViewer viewer)
1646: {
1647: PetscBool isascii;
1648: char str[50];
1649: const char *methodname[] = {
1650: "scaling & squaring, [m/m] Pade approximant (Higham)",
1651: "scaling & squaring, [6/6] Pade approximant",
1652: "scaling & squaring, subdiagonal Pade approximant (product form)",
1653: "scaling & squaring, subdiagonal Pade approximant (partial fraction)"
1654: };
1655: const int nmeth=PETSC_STATIC_ARRAY_LENGTH(methodname);
1657: PetscFunctionBegin;
1658: PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii));
1659: if (isascii) {
1660: if (fn->beta==(PetscScalar)1.0) {
1661: if (fn->alpha==(PetscScalar)1.0) PetscCall(PetscViewerASCIIPrintf(viewer," exponential: exp(x)\n"));
1662: else {
1663: PetscCall(SlepcSNPrintfScalar(str,sizeof(str),fn->alpha,PETSC_TRUE));
1664: PetscCall(PetscViewerASCIIPrintf(viewer," exponential: exp(%s*x)\n",str));
1665: }
1666: } else {
1667: PetscCall(SlepcSNPrintfScalar(str,sizeof(str),fn->beta,PETSC_TRUE));
1668: if (fn->alpha==(PetscScalar)1.0) PetscCall(PetscViewerASCIIPrintf(viewer," exponential: %s*exp(x)\n",str));
1669: else {
1670: PetscCall(PetscViewerASCIIPrintf(viewer," exponential: %s",str));
1671: PetscCall(PetscViewerASCIIUseTabs(viewer,PETSC_FALSE));
1672: PetscCall(SlepcSNPrintfScalar(str,sizeof(str),fn->alpha,PETSC_TRUE));
1673: PetscCall(PetscViewerASCIIPrintf(viewer,"*exp(%s*x)\n",str));
1674: PetscCall(PetscViewerASCIIUseTabs(viewer,PETSC_TRUE));
1675: }
1676: }
1677: if (fn->method<nmeth) PetscCall(PetscViewerASCIIPrintf(viewer," computing matrix functions with: %s\n",methodname[fn->method]));
1678: }
1679: PetscFunctionReturn(PETSC_SUCCESS);
1680: }
1682: /*MC
1683: FNEXP - FNEXP = "exp" - The exponential function $f(x)=e^x$.
1685: Level: beginner
1687: .seealso: [](sec:fn), `FN`, `FNType`, `FNSetType()`
1688: M*/
1690: SLEPC_EXTERN PetscErrorCode FNCreate_Exp(FN fn)
1691: {
1692: PetscFunctionBegin;
1693: fn->ops->evaluatefunction = FNEvaluateFunction_Exp;
1694: fn->ops->evaluatederivative = FNEvaluateDerivative_Exp;
1695: fn->ops->evaluatefunctionmat[0] = FNEvaluateFunctionMat_Exp_Higham;
1696: fn->ops->evaluatefunctionmat[1] = FNEvaluateFunctionMat_Exp_Pade;
1697: fn->ops->evaluatefunctionmat[2] = FNEvaluateFunctionMat_Exp_GuettelNakatsukasa; /* product form */
1698: fn->ops->evaluatefunctionmat[3] = FNEvaluateFunctionMat_Exp_GuettelNakatsukasa; /* partial fraction */
1699: #if PetscDefined(HAVE_CUDA)
1700: fn->ops->evaluatefunctionmatcuda[1] = FNEvaluateFunctionMat_Exp_Pade_CUDA;
1701: #if PetscDefined(HAVE_MAGMA)
1702: fn->ops->evaluatefunctionmatcuda[0] = FNEvaluateFunctionMat_Exp_Higham_CUDAm;
1703: fn->ops->evaluatefunctionmatcuda[1] = FNEvaluateFunctionMat_Exp_Pade_CUDAm;
1704: fn->ops->evaluatefunctionmatcuda[2] = FNEvaluateFunctionMat_Exp_GuettelNakatsukasa_CUDAm; /* product form */
1705: fn->ops->evaluatefunctionmatcuda[3] = FNEvaluateFunctionMat_Exp_GuettelNakatsukasa_CUDAm; /* partial fraction */
1706: #endif
1707: #endif
1708: fn->ops->view = FNView_Exp;
1709: PetscFunctionReturn(PETSC_SUCCESS);
1710: }