Actual source code: fninvsqrt.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: Inverse square root function x^(-1/2)
12: */
14: #include <slepc/private/fnimpl.h>
15: #include <slepcblaslapack.h>
17: static PetscErrorCode FNEvaluateFunction_Invsqrt(FN fn,PetscScalar x,PetscScalar *y)
18: {
19: PetscFunctionBegin;
20: PetscCheck(x!=0.0,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Function not defined in the requested value");
21: #if !PetscDefined(USE_COMPLEX)
22: PetscCheck(x>0.0,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Function not defined in the requested value");
23: #endif
24: *y = 1.0/PetscSqrtScalar(x);
25: PetscFunctionReturn(PETSC_SUCCESS);
26: }
28: static PetscErrorCode FNEvaluateDerivative_Invsqrt(FN fn,PetscScalar x,PetscScalar *y)
29: {
30: PetscFunctionBegin;
31: PetscCheck(x!=0.0,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Derivative not defined in the requested value");
32: #if !PetscDefined(USE_COMPLEX)
33: PetscCheck(x>0.0,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Derivative not defined in the requested value");
34: #endif
35: *y = -1.0/(2.0*PetscPowScalarReal(x,1.5));
36: PetscFunctionReturn(PETSC_SUCCESS);
37: }
39: static PetscErrorCode FNEvaluateFunctionMat_Invsqrt_Schur(FN fn,Mat A,Mat B)
40: {
41: PetscBLASInt n=0,ld,*ipiv;
42: PetscScalar *Ba,*Wa;
43: PetscInt m;
44: Mat W;
46: PetscFunctionBegin;
47: PetscCall(FN_AllocateWorkMat(fn,A,&W));
48: if (A!=B) PetscCall(MatCopy(A,B,SAME_NONZERO_PATTERN));
49: PetscCall(MatDenseGetArray(B,&Ba));
50: PetscCall(MatDenseGetArray(W,&Wa));
51: /* compute B = sqrtm(A) */
52: PetscCall(MatGetSize(A,&m,NULL));
53: PetscCall(PetscBLASIntCast(m,&n));
54: ld = n;
55: PetscCall(FNSqrtmSchur(fn,n,Ba,n,PETSC_FALSE));
56: /* compute B = A\B */
57: PetscCall(PetscMalloc1(ld,&ipiv));
58: PetscCallLAPACKInfo("LAPACKgesv",LAPACKgesv_(&n,&n,Wa,&ld,ipiv,Ba,&ld,&info));
59: PetscCall(PetscLogFlops(2.0*n*n*n/3.0+2.0*n*n*n));
60: PetscCall(PetscFree(ipiv));
61: PetscCall(MatDenseRestoreArray(W,&Wa));
62: PetscCall(MatDenseRestoreArray(B,&Ba));
63: PetscCall(FN_FreeWorkMat(fn,&W));
64: PetscFunctionReturn(PETSC_SUCCESS);
65: }
67: static PetscErrorCode FNEvaluateFunctionMatVec_Invsqrt_Schur(FN fn,Mat A,Vec v)
68: {
69: PetscBLASInt n=0,ld,*ipiv,one=1;
70: PetscScalar *Ba,*Wa;
71: PetscInt m;
72: Mat B,W;
74: PetscFunctionBegin;
75: PetscCall(FN_AllocateWorkMat(fn,A,&B));
76: PetscCall(FN_AllocateWorkMat(fn,A,&W));
77: PetscCall(MatDenseGetArray(B,&Ba));
78: PetscCall(MatDenseGetArray(W,&Wa));
79: /* compute B_1 = sqrtm(A)*e_1 */
80: PetscCall(MatGetSize(A,&m,NULL));
81: PetscCall(PetscBLASIntCast(m,&n));
82: ld = n;
83: PetscCall(FNSqrtmSchur(fn,n,Ba,n,PETSC_TRUE));
84: /* compute B_1 = A\B_1 */
85: PetscCall(PetscMalloc1(ld,&ipiv));
86: PetscCallLAPACKInfo("LAPACKgesv",LAPACKgesv_(&n,&one,Wa,&ld,ipiv,Ba,&ld,&info));
87: PetscCall(PetscFree(ipiv));
88: PetscCall(MatDenseRestoreArray(W,&Wa));
89: PetscCall(MatDenseRestoreArray(B,&Ba));
90: PetscCall(MatGetColumnVector(B,v,0));
91: PetscCall(FN_FreeWorkMat(fn,&W));
92: PetscCall(FN_FreeWorkMat(fn,&B));
93: PetscFunctionReturn(PETSC_SUCCESS);
94: }
96: static PetscErrorCode FNEvaluateFunctionMat_Invsqrt_DBP(FN fn,Mat A,Mat B)
97: {
98: PetscBLASInt n=0;
99: PetscScalar *T;
100: PetscInt m;
102: PetscFunctionBegin;
103: if (A!=B) PetscCall(MatCopy(A,B,SAME_NONZERO_PATTERN));
104: PetscCall(MatDenseGetArray(B,&T));
105: PetscCall(MatGetSize(A,&m,NULL));
106: PetscCall(PetscBLASIntCast(m,&n));
107: PetscCall(FNSqrtmDenmanBeavers(fn,n,T,n,PETSC_TRUE));
108: PetscCall(MatDenseRestoreArray(B,&T));
109: PetscFunctionReturn(PETSC_SUCCESS);
110: }
112: static PetscErrorCode FNEvaluateFunctionMat_Invsqrt_NS(FN fn,Mat A,Mat B)
113: {
114: PetscBLASInt n=0;
115: PetscScalar *T;
116: PetscInt m;
118: PetscFunctionBegin;
119: if (A!=B) PetscCall(MatCopy(A,B,SAME_NONZERO_PATTERN));
120: PetscCall(MatDenseGetArray(B,&T));
121: PetscCall(MatGetSize(A,&m,NULL));
122: PetscCall(PetscBLASIntCast(m,&n));
123: PetscCall(FNSqrtmNewtonSchulz(fn,n,T,n,PETSC_TRUE));
124: PetscCall(MatDenseRestoreArray(B,&T));
125: PetscFunctionReturn(PETSC_SUCCESS);
126: }
128: static PetscErrorCode FNEvaluateFunctionMat_Invsqrt_Sadeghi(FN fn,Mat A,Mat B)
129: {
130: PetscBLASInt n=0,ld,*ipiv;
131: PetscScalar *Ba,*Wa;
132: PetscInt m;
133: Mat W;
135: PetscFunctionBegin;
136: PetscCall(FN_AllocateWorkMat(fn,A,&W));
137: if (A!=B) PetscCall(MatCopy(A,B,SAME_NONZERO_PATTERN));
138: PetscCall(MatDenseGetArray(B,&Ba));
139: PetscCall(MatDenseGetArray(W,&Wa));
140: /* compute B = sqrtm(A) */
141: PetscCall(MatGetSize(A,&m,NULL));
142: PetscCall(PetscBLASIntCast(m,&n));
143: ld = n;
144: PetscCall(FNSqrtmSadeghi(fn,n,Ba,n));
145: /* compute B = A\B */
146: PetscCall(PetscMalloc1(ld,&ipiv));
147: PetscCallLAPACKInfo("LAPACKgesv",LAPACKgesv_(&n,&n,Wa,&ld,ipiv,Ba,&ld,&info));
148: PetscCall(PetscLogFlops(2.0*n*n*n/3.0+2.0*n*n*n));
149: PetscCall(PetscFree(ipiv));
150: PetscCall(MatDenseRestoreArray(W,&Wa));
151: PetscCall(MatDenseRestoreArray(B,&Ba));
152: PetscCall(FN_FreeWorkMat(fn,&W));
153: PetscFunctionReturn(PETSC_SUCCESS);
154: }
156: #if PetscDefined(HAVE_CUDA)
157: PetscErrorCode FNEvaluateFunctionMat_Invsqrt_NS_CUDA(FN fn,Mat A,Mat B)
158: {
159: PetscBLASInt n=0;
160: PetscScalar *Ba;
161: PetscInt m;
163: PetscFunctionBegin;
164: if (A!=B) PetscCall(MatCopy(A,B,SAME_NONZERO_PATTERN));
165: PetscCall(MatDenseCUDAGetArray(B,&Ba));
166: PetscCall(MatGetSize(A,&m,NULL));
167: PetscCall(PetscBLASIntCast(m,&n));
168: PetscCall(FNSqrtmNewtonSchulz_CUDA(fn,n,Ba,n,PETSC_TRUE));
169: PetscCall(MatDenseCUDARestoreArray(B,&Ba));
170: PetscFunctionReturn(PETSC_SUCCESS);
171: }
173: #if PetscDefined(HAVE_MAGMA)
174: #include <slepcmagma.h>
176: PetscErrorCode FNEvaluateFunctionMat_Invsqrt_DBP_CUDAm(FN fn,Mat A,Mat B)
177: {
178: PetscBLASInt n=0;
179: PetscScalar *T;
180: PetscInt m;
182: PetscFunctionBegin;
183: if (A!=B) PetscCall(MatCopy(A,B,SAME_NONZERO_PATTERN));
184: PetscCall(MatDenseCUDAGetArray(B,&T));
185: PetscCall(MatGetSize(A,&m,NULL));
186: PetscCall(PetscBLASIntCast(m,&n));
187: PetscCall(FNSqrtmDenmanBeavers_CUDAm(fn,n,T,n,PETSC_TRUE));
188: PetscCall(MatDenseCUDARestoreArray(B,&T));
189: PetscFunctionReturn(PETSC_SUCCESS);
190: }
192: PetscErrorCode FNEvaluateFunctionMat_Invsqrt_Sadeghi_CUDAm(FN fn,Mat A,Mat B)
193: {
194: PetscBLASInt n=0,ld,*ipiv;
195: PetscScalar *Ba,*Wa;
196: PetscInt m;
197: Mat W;
199: PetscFunctionBegin;
200: PetscCall(FN_AllocateWorkMat(fn,A,&W));
201: if (A!=B) PetscCall(MatCopy(A,B,SAME_NONZERO_PATTERN));
202: PetscCall(MatDenseCUDAGetArray(B,&Ba));
203: PetscCall(MatDenseCUDAGetArray(W,&Wa));
204: /* compute B = sqrtm(A) */
205: PetscCall(MatGetSize(A,&m,NULL));
206: PetscCall(PetscBLASIntCast(m,&n));
207: ld = n;
208: PetscCall(FNSqrtmSadeghi_CUDAm(fn,n,Ba,n));
209: /* compute B = A\B */
210: PetscCall(SlepcMagmaInit());
211: PetscCall(PetscMalloc1(ld,&ipiv));
212: PetscCallMAGMA(magma_xgesv_gpu,n,n,Wa,ld,ipiv,Ba,ld);
213: PetscCall(PetscLogFlops(2.0*n*n*n/3.0+2.0*n*n*n));
214: PetscCall(PetscFree(ipiv));
215: PetscCall(MatDenseCUDARestoreArray(W,&Wa));
216: PetscCall(MatDenseCUDARestoreArray(B,&Ba));
217: PetscCall(FN_FreeWorkMat(fn,&W));
218: PetscFunctionReturn(PETSC_SUCCESS);
219: }
220: #endif /* PETSC_HAVE_MAGMA */
221: #endif /* PETSC_HAVE_CUDA */
223: static PetscErrorCode FNView_Invsqrt(FN fn,PetscViewer viewer)
224: {
225: PetscBool isascii;
226: char str[50];
227: const char *methodname[] = {
228: "Schur method for inv(A)*sqrtm(A)",
229: "Denman-Beavers (product form)",
230: "Newton-Schulz iteration",
231: "Sadeghi iteration"
232: };
233: const int nmeth=PETSC_STATIC_ARRAY_LENGTH(methodname);
235: PetscFunctionBegin;
236: PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii));
237: if (isascii) {
238: if (fn->beta==(PetscScalar)1.0) {
239: if (fn->alpha==(PetscScalar)1.0) PetscCall(PetscViewerASCIIPrintf(viewer," inverse square root: x^(-1/2)\n"));
240: else {
241: PetscCall(SlepcSNPrintfScalar(str,sizeof(str),fn->alpha,PETSC_TRUE));
242: PetscCall(PetscViewerASCIIPrintf(viewer," inverse square root: (%s*x)^(-1/2)\n",str));
243: }
244: } else {
245: PetscCall(SlepcSNPrintfScalar(str,sizeof(str),fn->beta,PETSC_TRUE));
246: if (fn->alpha==(PetscScalar)1.0) PetscCall(PetscViewerASCIIPrintf(viewer," inverse square root: %s*x^(-1/2)\n",str));
247: else {
248: PetscCall(PetscViewerASCIIPrintf(viewer," inverse square root: %s",str));
249: PetscCall(PetscViewerASCIIUseTabs(viewer,PETSC_FALSE));
250: PetscCall(SlepcSNPrintfScalar(str,sizeof(str),fn->alpha,PETSC_TRUE));
251: PetscCall(PetscViewerASCIIPrintf(viewer,"*(%s*x)^(-1/2)\n",str));
252: PetscCall(PetscViewerASCIIUseTabs(viewer,PETSC_TRUE));
253: }
254: }
255: if (fn->method<nmeth) PetscCall(PetscViewerASCIIPrintf(viewer," computing matrix functions with: %s\n",methodname[fn->method]));
256: }
257: PetscFunctionReturn(PETSC_SUCCESS);
258: }
260: /*MC
261: FNINVSQRT - FNINVSQRT = "invsqrt" - The inverse square root function $f(x)=x^{-\frac{1}{2}}$.
263: Level: beginner
265: .seealso: [](sec:fn), `FN`, `FNType`, `FNSetType()`
266: M*/
268: SLEPC_EXTERN PetscErrorCode FNCreate_Invsqrt(FN fn)
269: {
270: PetscFunctionBegin;
271: fn->ops->evaluatefunction = FNEvaluateFunction_Invsqrt;
272: fn->ops->evaluatederivative = FNEvaluateDerivative_Invsqrt;
273: fn->ops->evaluatefunctionmat[0] = FNEvaluateFunctionMat_Invsqrt_Schur;
274: fn->ops->evaluatefunctionmat[1] = FNEvaluateFunctionMat_Invsqrt_DBP;
275: fn->ops->evaluatefunctionmat[2] = FNEvaluateFunctionMat_Invsqrt_NS;
276: fn->ops->evaluatefunctionmat[3] = FNEvaluateFunctionMat_Invsqrt_Sadeghi;
277: #if PetscDefined(HAVE_CUDA)
278: fn->ops->evaluatefunctionmatcuda[2] = FNEvaluateFunctionMat_Invsqrt_NS_CUDA;
279: #if PetscDefined(HAVE_MAGMA)
280: fn->ops->evaluatefunctionmatcuda[1] = FNEvaluateFunctionMat_Invsqrt_DBP_CUDAm;
281: fn->ops->evaluatefunctionmatcuda[3] = FNEvaluateFunctionMat_Invsqrt_Sadeghi_CUDAm;
282: #endif /* PETSC_HAVE_MAGMA */
283: #endif /* PETSC_HAVE_CUDA */
284: fn->ops->evaluatefunctionmatvec[0] = FNEvaluateFunctionMatVec_Invsqrt_Schur;
285: fn->ops->view = FNView_Invsqrt;
286: PetscFunctionReturn(PETSC_SUCCESS);
287: }