Actual source code: fninvsqrt.c
slepc-3.11.2 2019-07-30
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-2019, 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> /*I "slepcfn.h" I*/
15: #include <slepcblaslapack.h>
17: PetscErrorCode FNEvaluateFunction_Invsqrt(FN fn,PetscScalar x,PetscScalar *y)
18: {
20: if (x==0.0) SETERRQ(PETSC_COMM_SELF,1,"Function not defined in the requested value");
21: #if !defined(PETSC_USE_COMPLEX)
22: if (x<0.0) SETERRQ(PETSC_COMM_SELF,1,"Function not defined in the requested value");
23: #endif
24: *y = 1.0/PetscSqrtScalar(x);
25: return(0);
26: }
28: PetscErrorCode FNEvaluateDerivative_Invsqrt(FN fn,PetscScalar x,PetscScalar *y)
29: {
31: if (x==0.0) SETERRQ(PETSC_COMM_SELF,1,"Derivative not defined in the requested value");
32: #if !defined(PETSC_USE_COMPLEX)
33: if (x<0.0) SETERRQ(PETSC_COMM_SELF,1,"Derivative not defined in the requested value");
34: #endif
35: *y = -1.0/(2.0*PetscPowScalarReal(x,1.5));
36: return(0);
37: }
39: PetscErrorCode FNEvaluateFunctionMat_Invsqrt_Schur(FN fn,Mat A,Mat B)
40: {
41: #if defined(PETSC_MISSING_LAPACK_GESV)
43: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"GESV - Lapack routine is unavailable");
44: #else
46: PetscBLASInt n,ld,*ipiv,info;
47: PetscScalar *Ba,*Wa;
48: PetscInt m;
49: Mat W;
52: FN_AllocateWorkMat(fn,A,&W);
53: if (A!=B) { MatCopy(A,B,SAME_NONZERO_PATTERN); }
54: MatDenseGetArray(B,&Ba);
55: MatDenseGetArray(W,&Wa);
56: /* compute B = sqrtm(A) */
57: MatGetSize(A,&m,NULL);
58: PetscBLASIntCast(m,&n);
59: ld = n;
60: SlepcSqrtmSchur(n,Ba,n,PETSC_FALSE);
61: /* compute B = A\B */
62: PetscMalloc1(ld,&ipiv);
63: PetscStackCallBLAS("LAPACKgesv",LAPACKgesv_(&n,&n,Wa,&ld,ipiv,Ba,&ld,&info));
64: SlepcCheckLapackInfo("gesv",info);
65: PetscLogFlops(2.0*n*n*n/3.0+2.0*n*n*n);
66: PetscFree(ipiv);
67: MatDenseRestoreArray(W,&Wa);
68: MatDenseRestoreArray(B,&Ba);
69: FN_FreeWorkMat(fn,&W);
70: return(0);
71: #endif
72: }
74: PetscErrorCode FNEvaluateFunctionMatVec_Invsqrt_Schur(FN fn,Mat A,Vec v)
75: {
76: #if defined(PETSC_MISSING_LAPACK_GESV)
78: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"GESV - Lapack routine is unavailable");
79: #else
81: PetscBLASInt n,ld,*ipiv,info,one=1;
82: PetscScalar *Ba,*Wa;
83: PetscInt m;
84: Mat B,W;
87: FN_AllocateWorkMat(fn,A,&B);
88: FN_AllocateWorkMat(fn,A,&W);
89: MatDenseGetArray(B,&Ba);
90: MatDenseGetArray(W,&Wa);
91: /* compute B_1 = sqrtm(A)*e_1 */
92: MatGetSize(A,&m,NULL);
93: PetscBLASIntCast(m,&n);
94: ld = n;
95: SlepcSqrtmSchur(n,Ba,n,PETSC_TRUE);
96: /* compute B_1 = A\B_1 */
97: PetscMalloc1(ld,&ipiv);
98: PetscStackCallBLAS("LAPACKgesv",LAPACKgesv_(&n,&one,Wa,&ld,ipiv,Ba,&ld,&info));
99: SlepcCheckLapackInfo("gesv",info);
100: PetscFree(ipiv);
101: MatDenseRestoreArray(W,&Wa);
102: MatDenseRestoreArray(B,&Ba);
103: MatGetColumnVector(B,v,0);
104: FN_FreeWorkMat(fn,&W);
105: FN_FreeWorkMat(fn,&B);
106: return(0);
107: #endif
108: }
110: PetscErrorCode FNEvaluateFunctionMat_Invsqrt_DBP(FN fn,Mat A,Mat B)
111: {
113: PetscBLASInt n;
114: PetscScalar *T;
115: PetscInt m;
118: if (A!=B) { MatCopy(A,B,SAME_NONZERO_PATTERN); }
119: MatDenseGetArray(B,&T);
120: MatGetSize(A,&m,NULL);
121: PetscBLASIntCast(m,&n);
122: SlepcSqrtmDenmanBeavers(n,T,n,PETSC_TRUE);
123: MatDenseRestoreArray(B,&T);
124: return(0);
125: }
127: PetscErrorCode FNEvaluateFunctionMat_Invsqrt_NS(FN fn,Mat A,Mat B)
128: {
130: PetscBLASInt n;
131: PetscScalar *T;
132: PetscInt m;
135: if (A!=B) { MatCopy(A,B,SAME_NONZERO_PATTERN); }
136: MatDenseGetArray(B,&T);
137: MatGetSize(A,&m,NULL);
138: PetscBLASIntCast(m,&n);
139: SlepcSqrtmNewtonSchulz(n,T,n,PETSC_TRUE);
140: MatDenseRestoreArray(B,&T);
141: return(0);
142: }
144: PetscErrorCode FNView_Invsqrt(FN fn,PetscViewer viewer)
145: {
147: PetscBool isascii;
148: char str[50];
149: const char *methodname[] = {
150: "Schur method for inv(A)*sqrtm(A)",
151: "Denman-Beavers (product form)",
152: "Newton-Schulz iteration"
153: };
154: const int nmeth=sizeof(methodname)/sizeof(methodname[0]);
157: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
158: if (isascii) {
159: if (fn->beta==(PetscScalar)1.0) {
160: if (fn->alpha==(PetscScalar)1.0) {
161: PetscViewerASCIIPrintf(viewer," Inverse square root: x^(-1/2)\n");
162: } else {
163: SlepcSNPrintfScalar(str,50,fn->alpha,PETSC_TRUE);
164: PetscViewerASCIIPrintf(viewer," Inverse square root: (%s*x)^(-1/2)\n",str);
165: }
166: } else {
167: SlepcSNPrintfScalar(str,50,fn->beta,PETSC_TRUE);
168: if (fn->alpha==(PetscScalar)1.0) {
169: PetscViewerASCIIPrintf(viewer," Inverse square root: %s*x^(-1/2)\n",str);
170: } else {
171: PetscViewerASCIIPrintf(viewer," Inverse square root: %s",str);
172: PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
173: SlepcSNPrintfScalar(str,50,fn->alpha,PETSC_TRUE);
174: PetscViewerASCIIPrintf(viewer,"*(%s*x)^(-1/2)\n",str);
175: PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
176: }
177: }
178: if (fn->method<nmeth) {
179: PetscViewerASCIIPrintf(viewer," computing matrix functions with: %s\n",methodname[fn->method]);
180: }
181: }
182: return(0);
183: }
185: SLEPC_EXTERN PetscErrorCode FNCreate_Invsqrt(FN fn)
186: {
188: fn->ops->evaluatefunction = FNEvaluateFunction_Invsqrt;
189: fn->ops->evaluatederivative = FNEvaluateDerivative_Invsqrt;
190: fn->ops->evaluatefunctionmat[0] = FNEvaluateFunctionMat_Invsqrt_Schur;
191: fn->ops->evaluatefunctionmat[1] = FNEvaluateFunctionMat_Invsqrt_DBP;
192: fn->ops->evaluatefunctionmat[2] = FNEvaluateFunctionMat_Invsqrt_NS;
193: fn->ops->evaluatefunctionmatvec[0] = FNEvaluateFunctionMatVec_Invsqrt_Schur;
194: fn->ops->view = FNView_Invsqrt;
195: return(0);
196: }