Actual source code: fninvsqrt.c

slepc-3.11.2 2019-07-30
Report Typos and Errors
  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: }