Actual source code: fninvsqrt.c
slepc-3.14.1 2020-12-08
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-2020, 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: 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: {
42: PetscBLASInt n,ld,*ipiv,info;
43: PetscScalar *Ba,*Wa;
44: PetscInt m;
45: Mat W;
48: FN_AllocateWorkMat(fn,A,&W);
49: if (A!=B) { MatCopy(A,B,SAME_NONZERO_PATTERN); }
50: MatDenseGetArray(B,&Ba);
51: MatDenseGetArray(W,&Wa);
52: /* compute B = sqrtm(A) */
53: MatGetSize(A,&m,NULL);
54: PetscBLASIntCast(m,&n);
55: ld = n;
56: SlepcSqrtmSchur(n,Ba,n,PETSC_FALSE);
57: /* compute B = A\B */
58: PetscMalloc1(ld,&ipiv);
59: PetscStackCallBLAS("LAPACKgesv",LAPACKgesv_(&n,&n,Wa,&ld,ipiv,Ba,&ld,&info));
60: SlepcCheckLapackInfo("gesv",info);
61: PetscLogFlops(2.0*n*n*n/3.0+2.0*n*n*n);
62: PetscFree(ipiv);
63: MatDenseRestoreArray(W,&Wa);
64: MatDenseRestoreArray(B,&Ba);
65: FN_FreeWorkMat(fn,&W);
66: return(0);
67: }
69: PetscErrorCode FNEvaluateFunctionMatVec_Invsqrt_Schur(FN fn,Mat A,Vec v)
70: {
72: PetscBLASInt n,ld,*ipiv,info,one=1;
73: PetscScalar *Ba,*Wa;
74: PetscInt m;
75: Mat B,W;
78: FN_AllocateWorkMat(fn,A,&B);
79: FN_AllocateWorkMat(fn,A,&W);
80: MatDenseGetArray(B,&Ba);
81: MatDenseGetArray(W,&Wa);
82: /* compute B_1 = sqrtm(A)*e_1 */
83: MatGetSize(A,&m,NULL);
84: PetscBLASIntCast(m,&n);
85: ld = n;
86: SlepcSqrtmSchur(n,Ba,n,PETSC_TRUE);
87: /* compute B_1 = A\B_1 */
88: PetscMalloc1(ld,&ipiv);
89: PetscStackCallBLAS("LAPACKgesv",LAPACKgesv_(&n,&one,Wa,&ld,ipiv,Ba,&ld,&info));
90: SlepcCheckLapackInfo("gesv",info);
91: PetscFree(ipiv);
92: MatDenseRestoreArray(W,&Wa);
93: MatDenseRestoreArray(B,&Ba);
94: MatGetColumnVector(B,v,0);
95: FN_FreeWorkMat(fn,&W);
96: FN_FreeWorkMat(fn,&B);
97: return(0);
98: }
100: PetscErrorCode FNEvaluateFunctionMat_Invsqrt_DBP(FN fn,Mat A,Mat B)
101: {
103: PetscBLASInt n = 0;
104: PetscScalar *T;
105: PetscInt m;
108: if (A!=B) { MatCopy(A,B,SAME_NONZERO_PATTERN); }
109: MatDenseGetArray(B,&T);
110: MatGetSize(A,&m,NULL);
111: PetscBLASIntCast(m,&n);
112: SlepcSqrtmDenmanBeavers(n,T,n,PETSC_TRUE);
113: MatDenseRestoreArray(B,&T);
114: return(0);
115: }
117: PetscErrorCode FNEvaluateFunctionMat_Invsqrt_NS(FN fn,Mat A,Mat B)
118: {
120: PetscBLASInt n = 0;
121: PetscScalar *T;
122: PetscInt m;
125: if (A!=B) { MatCopy(A,B,SAME_NONZERO_PATTERN); }
126: MatDenseGetArray(B,&T);
127: MatGetSize(A,&m,NULL);
128: PetscBLASIntCast(m,&n);
129: SlepcSqrtmNewtonSchulz(n,T,n,PETSC_TRUE);
130: MatDenseRestoreArray(B,&T);
131: return(0);
132: }
134: PetscErrorCode FNView_Invsqrt(FN fn,PetscViewer viewer)
135: {
137: PetscBool isascii;
138: char str[50];
139: const char *methodname[] = {
140: "Schur method for inv(A)*sqrtm(A)",
141: "Denman-Beavers (product form)",
142: "Newton-Schulz iteration"
143: };
144: const int nmeth=sizeof(methodname)/sizeof(methodname[0]);
147: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
148: if (isascii) {
149: if (fn->beta==(PetscScalar)1.0) {
150: if (fn->alpha==(PetscScalar)1.0) {
151: PetscViewerASCIIPrintf(viewer," Inverse square root: x^(-1/2)\n");
152: } else {
153: SlepcSNPrintfScalar(str,sizeof(str),fn->alpha,PETSC_TRUE);
154: PetscViewerASCIIPrintf(viewer," Inverse square root: (%s*x)^(-1/2)\n",str);
155: }
156: } else {
157: SlepcSNPrintfScalar(str,sizeof(str),fn->beta,PETSC_TRUE);
158: if (fn->alpha==(PetscScalar)1.0) {
159: PetscViewerASCIIPrintf(viewer," Inverse square root: %s*x^(-1/2)\n",str);
160: } else {
161: PetscViewerASCIIPrintf(viewer," Inverse square root: %s",str);
162: PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
163: SlepcSNPrintfScalar(str,sizeof(str),fn->alpha,PETSC_TRUE);
164: PetscViewerASCIIPrintf(viewer,"*(%s*x)^(-1/2)\n",str);
165: PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
166: }
167: }
168: if (fn->method<nmeth) {
169: PetscViewerASCIIPrintf(viewer," computing matrix functions with: %s\n",methodname[fn->method]);
170: }
171: }
172: return(0);
173: }
175: SLEPC_EXTERN PetscErrorCode FNCreate_Invsqrt(FN fn)
176: {
178: fn->ops->evaluatefunction = FNEvaluateFunction_Invsqrt;
179: fn->ops->evaluatederivative = FNEvaluateDerivative_Invsqrt;
180: fn->ops->evaluatefunctionmat[0] = FNEvaluateFunctionMat_Invsqrt_Schur;
181: fn->ops->evaluatefunctionmat[1] = FNEvaluateFunctionMat_Invsqrt_DBP;
182: fn->ops->evaluatefunctionmat[2] = FNEvaluateFunctionMat_Invsqrt_NS;
183: fn->ops->evaluatefunctionmatvec[0] = FNEvaluateFunctionMatVec_Invsqrt_Schur;
184: fn->ops->view = FNView_Invsqrt;
185: return(0);
186: }