Actual source code: fnsqrt.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:    Square root function  sqrt(x)
 12: */

 14: #include <slepc/private/fnimpl.h>
 15: #include <slepcblaslapack.h>

 17: static PetscErrorCode FNEvaluateFunction_Sqrt(FN fn,PetscScalar x,PetscScalar *y)
 18: {
 19:   PetscFunctionBegin;
 20: #if !PetscDefined(USE_COMPLEX)
 21:   PetscCheck(x>=0.0,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Function not defined in the requested value");
 22: #endif
 23:   *y = PetscSqrtScalar(x);
 24:   PetscFunctionReturn(PETSC_SUCCESS);
 25: }

 27: static PetscErrorCode FNEvaluateDerivative_Sqrt(FN fn,PetscScalar x,PetscScalar *y)
 28: {
 29:   PetscFunctionBegin;
 30:   PetscCheck(x!=0.0,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Derivative not defined in the requested value");
 31: #if !PetscDefined(USE_COMPLEX)
 32:   PetscCheck(x>0.0,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Derivative not defined in the requested value");
 33: #endif
 34:   *y = 1.0/(2.0*PetscSqrtScalar(x));
 35:   PetscFunctionReturn(PETSC_SUCCESS);
 36: }

 38: static PetscErrorCode FNEvaluateFunctionMat_Sqrt_Schur(FN fn,Mat A,Mat B)
 39: {
 40:   PetscBLASInt   n=0;
 41:   PetscScalar    *T;
 42:   PetscInt       m;

 44:   PetscFunctionBegin;
 45:   if (A!=B) PetscCall(MatCopy(A,B,SAME_NONZERO_PATTERN));
 46:   PetscCall(MatDenseGetArray(B,&T));
 47:   PetscCall(MatGetSize(A,&m,NULL));
 48:   PetscCall(PetscBLASIntCast(m,&n));
 49:   PetscCall(FNSqrtmSchur(fn,n,T,n,PETSC_FALSE));
 50:   PetscCall(MatDenseRestoreArray(B,&T));
 51:   PetscFunctionReturn(PETSC_SUCCESS);
 52: }

 54: static PetscErrorCode FNEvaluateFunctionMatVec_Sqrt_Schur(FN fn,Mat A,Vec v)
 55: {
 56:   PetscBLASInt   n=0;
 57:   PetscScalar    *T;
 58:   PetscInt       m;
 59:   Mat            B;

 61:   PetscFunctionBegin;
 62:   PetscCall(FN_AllocateWorkMat(fn,A,&B));
 63:   PetscCall(MatDenseGetArray(B,&T));
 64:   PetscCall(MatGetSize(A,&m,NULL));
 65:   PetscCall(PetscBLASIntCast(m,&n));
 66:   PetscCall(FNSqrtmSchur(fn,n,T,n,PETSC_TRUE));
 67:   PetscCall(MatDenseRestoreArray(B,&T));
 68:   PetscCall(MatGetColumnVector(B,v,0));
 69:   PetscCall(FN_FreeWorkMat(fn,&B));
 70:   PetscFunctionReturn(PETSC_SUCCESS);
 71: }

 73: static PetscErrorCode FNEvaluateFunctionMat_Sqrt_DBP(FN fn,Mat A,Mat B)
 74: {
 75:   PetscBLASInt   n=0;
 76:   PetscScalar    *T;
 77:   PetscInt       m;

 79:   PetscFunctionBegin;
 80:   if (A!=B) PetscCall(MatCopy(A,B,SAME_NONZERO_PATTERN));
 81:   PetscCall(MatDenseGetArray(B,&T));
 82:   PetscCall(MatGetSize(A,&m,NULL));
 83:   PetscCall(PetscBLASIntCast(m,&n));
 84:   PetscCall(FNSqrtmDenmanBeavers(fn,n,T,n,PETSC_FALSE));
 85:   PetscCall(MatDenseRestoreArray(B,&T));
 86:   PetscFunctionReturn(PETSC_SUCCESS);
 87: }

 89: static PetscErrorCode FNEvaluateFunctionMat_Sqrt_NS(FN fn,Mat A,Mat B)
 90: {
 91:   PetscBLASInt   n=0;
 92:   PetscScalar    *Ba;
 93:   PetscInt       m;

 95:   PetscFunctionBegin;
 96:   if (A!=B) PetscCall(MatCopy(A,B,SAME_NONZERO_PATTERN));
 97:   PetscCall(MatDenseGetArray(B,&Ba));
 98:   PetscCall(MatGetSize(A,&m,NULL));
 99:   PetscCall(PetscBLASIntCast(m,&n));
100:   PetscCall(FNSqrtmNewtonSchulz(fn,n,Ba,n,PETSC_FALSE));
101:   PetscCall(MatDenseRestoreArray(B,&Ba));
102:   PetscFunctionReturn(PETSC_SUCCESS);
103: }

105: #define MAXIT 50

107: /*
108:    Computes the principal square root of the matrix A using the
109:    Sadeghi iteration. A is overwritten with sqrtm(A).
110:  */
111: PetscErrorCode FNSqrtmSadeghi(FN fn,PetscBLASInt n,PetscScalar *A,PetscBLASInt ld)
112: {
113:   PetscScalar    *M,*M2,*G,*X=A,*work,work1,sqrtnrm;
114:   PetscScalar    szero=0.0,sone=1.0,smfive=-5.0,s1d16=1.0/16.0;
115:   PetscReal      tol,Mres=0.0,nrm,rwork[1],done=1.0;
116:   PetscInt       i,it;
117:   PetscBLASInt   N,*piv=NULL,lwork=0,query=-1,one=1,zero=0;
118:   PetscBool      converged=PETSC_FALSE;
119:   unsigned int   ftz;

121:   PetscFunctionBegin;
122:   N = n*n;
123:   tol = PetscSqrtReal((PetscReal)n)*PETSC_MACHINE_EPSILON/2;
124:   PetscCall(SlepcSetFlushToZero(&ftz));

126:   /* query work size */
127:   PetscCallLAPACKInfo("LAPACKgetri",LAPACKgetri_(&n,A,&ld,piv,&work1,&query,&info));
128:   PetscCall(PetscBLASIntCast((PetscInt)PetscRealPart(work1),&lwork));

130:   PetscCall(PetscMalloc5(N,&M,N,&M2,N,&G,lwork,&work,n,&piv));
131:   PetscCall(PetscArraycpy(M,A,N));

133:   /* scale M */
134:   nrm = LAPACKlange_("fro",&n,&n,M,&n,rwork);
135:   if (nrm>1.0) {
136:     sqrtnrm = PetscSqrtReal(nrm);
137:     PetscCallLAPACKInfo("LAPACKlascl",LAPACKlascl_("G",&zero,&zero,&nrm,&done,&N,&one,M,&N,&info));
138:     tol *= nrm;
139:   }
140:   PetscCall(PetscInfo(fn,"||A||_F = %g, new tol: %g\n",(double)nrm,(double)tol));

142:   /* X = I */
143:   PetscCall(PetscArrayzero(X,N));
144:   for (i=0;i<n;i++) X[i+i*ld] = 1.0;

146:   for (it=0;it<MAXIT && !converged;it++) {

148:     /* G = (5/16)*I + (1/16)*M*(15*I-5*M+M*M) */
149:     PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&n,&n,&n,&sone,M,&ld,M,&ld,&szero,M2,&ld));
150:     PetscCallBLAS("BLASaxpy",BLASaxpy_(&N,&smfive,M,&one,M2,&one));
151:     for (i=0;i<n;i++) M2[i+i*ld] += 15.0;
152:     PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&n,&n,&n,&s1d16,M,&ld,M2,&ld,&szero,G,&ld));
153:     for (i=0;i<n;i++) G[i+i*ld] += 5.0/16.0;

155:     /* X = X*G */
156:     PetscCall(PetscArraycpy(M2,X,N));
157:     PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&n,&n,&n,&sone,M2,&ld,G,&ld,&szero,X,&ld));

159:     /* M = M*inv(G*G) */
160:     PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&n,&n,&n,&sone,G,&ld,G,&ld,&szero,M2,&ld));
161:     PetscCallLAPACKInfo("LAPACKgetrf",LAPACKgetrf_(&n,&n,M2,&ld,piv,&info));
162:     PetscCallLAPACKInfo("LAPACKgetri",LAPACKgetri_(&n,M2,&ld,piv,work,&lwork,&info));

164:     PetscCall(PetscArraycpy(G,M,N));
165:     PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&n,&n,&n,&sone,G,&ld,M2,&ld,&szero,M,&ld));

167:     /* check ||I-M|| */
168:     PetscCall(PetscArraycpy(M2,M,N));
169:     for (i=0;i<n;i++) M2[i+i*ld] -= 1.0;
170:     Mres = LAPACKlange_("fro",&n,&n,M2,&n,rwork);
171:     PetscCheck(!PetscIsNanReal(Mres),PETSC_COMM_SELF,PETSC_ERR_FP,"The computed norm is not-a-number");
172:     if (Mres<=tol) converged = PETSC_TRUE;
173:     PetscCall(PetscInfo(fn,"it: %" PetscInt_FMT " res: %g\n",it,(double)Mres));
174:     PetscCall(PetscLogFlops(8.0*n*n*n+2.0*n*n+2.0*n*n*n/3.0+4.0*n*n*n/3.0+2.0*n*n*n+2.0*n*n));
175:   }

177:   PetscCheck(Mres<=tol,PETSC_COMM_SELF,PETSC_ERR_LIB,"SQRTM not converged after %d iterations",MAXIT);

179:   /* undo scaling */
180:   if (nrm>1.0) PetscCallBLAS("BLASscal",BLASscal_(&N,&sqrtnrm,A,&one));

182:   PetscCall(PetscFree5(M,M2,G,work,piv));
183:   PetscCall(SlepcResetFlushToZero(&ftz));
184:   PetscFunctionReturn(PETSC_SUCCESS);
185: }

187: #if PetscDefined(HAVE_CUDA)
188: #include "../src/sys/classes/fn/impls/cuda/fnutilcuda.h"
189: #include <slepccupmblas.h>

191: #if PetscDefined(HAVE_MAGMA)
192: #include <slepcmagma.h>

194: /*
195:  * Matrix square root by Sadeghi iteration. CUDA version.
196:  * Computes the principal square root of the matrix A using the
197:  * Sadeghi iteration. A is overwritten with sqrtm(A).
198:  */
199: PetscErrorCode FNSqrtmSadeghi_CUDAm(FN fn,PetscBLASInt n,PetscScalar *d_A,PetscBLASInt ld)
200: {
201:   PetscScalar        *d_M,*d_M2,*d_G,*d_work,alpha;
202:   const PetscScalar  szero=0.0,sone=1.0,smfive=-5.0,s15=15.0,s1d16=1.0/16.0;
203:   PetscReal          tol,Mres=0.0,nrm,sqrtnrm=1.0;
204:   PetscInt           it,nb,lwork;
205:   PetscBLASInt       *piv,N;
206:   const PetscBLASInt one=1;
207:   PetscBool          converged=PETSC_FALSE;
208:   cublasHandle_t     cublasv2handle;

210:   PetscFunctionBegin;
211:   PetscCall(PetscDeviceInitialize(PETSC_DEVICE_CUDA)); /* For CUDA event timers */
212:   PetscCall(PetscCUBLASGetHandle(&cublasv2handle));
213:   PetscCall(SlepcMagmaInit());
214:   N = n*n;
215:   tol = PetscSqrtReal((PetscReal)n)*PETSC_MACHINE_EPSILON/2;

217:   PetscCall(PetscMalloc1(n,&piv));
218:   PetscCallCUDA(cudaMalloc((void **)&d_M,sizeof(PetscScalar)*N));
219:   PetscCallCUDA(cudaMalloc((void **)&d_M2,sizeof(PetscScalar)*N));
220:   PetscCallCUDA(cudaMalloc((void **)&d_G,sizeof(PetscScalar)*N));

222:   nb = magma_get_xgetri_nb(n);
223:   lwork = nb*n;
224:   PetscCallCUDA(cudaMalloc((void **)&d_work,sizeof(PetscScalar)*lwork));
225:   PetscCall(PetscLogGpuTimeBegin());

227:   /* M = A */
228:   PetscCallCUDA(cudaMemcpy(d_M,d_A,sizeof(PetscScalar)*N,cudaMemcpyDeviceToDevice));

230:   /* scale M */
231:   PetscCallCUBLAS(cublasXnrm2(cublasv2handle,N,d_M,one,&nrm));
232:   if (nrm>1.0) {
233:     sqrtnrm = PetscSqrtReal(nrm);
234:     alpha = 1.0/nrm;
235:     PetscCallCUBLAS(cublasXscal(cublasv2handle,N,&alpha,d_M,one));
236:     tol *= nrm;
237:   }
238:   PetscCall(PetscInfo(fn,"||A||_F = %g, new tol: %g\n",(double)nrm,(double)tol));

240:   /* X = I */
241:   PetscCallCUDA(cudaMemset(d_A,0,sizeof(PetscScalar)*N));
242:   PetscCall(set_diagonal(n,d_A,ld,sone));

244:   for (it=0;it<MAXIT && !converged;it++) {

246:     /* G = (5/16)*I + (1/16)*M*(15*I-5*M+M*M) */
247:     PetscCallCUBLAS(cublasXgemm(cublasv2handle,CUBLAS_OP_N,CUBLAS_OP_N,n,n,n,&sone,d_M,ld,d_M,ld,&szero,d_M2,ld));
248:     PetscCallCUBLAS(cublasXaxpy(cublasv2handle,N,&smfive,d_M,one,d_M2,one));
249:     PetscCall(shift_diagonal(n,d_M2,ld,s15));
250:     PetscCallCUBLAS(cublasXgemm(cublasv2handle,CUBLAS_OP_N,CUBLAS_OP_N,n,n,n,&s1d16,d_M,ld,d_M2,ld,&szero,d_G,ld));
251:     PetscCall(shift_diagonal(n,d_G,ld,5.0/16.0));

253:     /* X = X*G */
254:     PetscCallCUDA(cudaMemcpy(d_M2,d_A,sizeof(PetscScalar)*N,cudaMemcpyDeviceToDevice));
255:     PetscCallCUBLAS(cublasXgemm(cublasv2handle,CUBLAS_OP_N,CUBLAS_OP_N,n,n,n,&sone,d_M2,ld,d_G,ld,&szero,d_A,ld));

257:     /* M = M*inv(G*G) */
258:     PetscCallCUBLAS(cublasXgemm(cublasv2handle,CUBLAS_OP_N,CUBLAS_OP_N,n,n,n,&sone,d_G,ld,d_G,ld,&szero,d_M2,ld));
259:     /* magma */
260:     PetscCallMAGMA(magma_xgetrf_gpu,n,n,d_M2,ld,piv);
261:     PetscCallMAGMA(magma_xgetri_gpu,n,d_M2,ld,piv,d_work,lwork);
262:     /* magma */
263:     PetscCallCUDA(cudaMemcpy(d_G,d_M,sizeof(PetscScalar)*N,cudaMemcpyDeviceToDevice));
264:     PetscCallCUBLAS(cublasXgemm(cublasv2handle,CUBLAS_OP_N,CUBLAS_OP_N,n,n,n,&sone,d_G,ld,d_M2,ld,&szero,d_M,ld));

266:     /* check ||I-M|| */
267:     PetscCallCUDA(cudaMemcpy(d_M2,d_M,sizeof(PetscScalar)*N,cudaMemcpyDeviceToDevice));
268:     PetscCall(shift_diagonal(n,d_M2,ld,-1.0));
269:     PetscCallCUBLAS(cublasXnrm2(cublasv2handle,N,d_M2,one,&Mres));
270:     PetscCheck(!PetscIsNanReal(Mres),PETSC_COMM_SELF,PETSC_ERR_FP,"The computed norm is not-a-number");
271:     if (Mres<=tol) converged = PETSC_TRUE;
272:     PetscCall(PetscInfo(fn,"it: %" PetscInt_FMT " res: %g\n",it,(double)Mres));
273:     PetscCall(PetscLogGpuFlops(8.0*n*n*n+2.0*n*n+2.0*n*n*n/3.0+4.0*n*n*n/3.0+2.0*n*n*n+2.0*n*n));
274:   }

276:   PetscCheck(Mres<=tol,PETSC_COMM_SELF,PETSC_ERR_LIB,"SQRTM not converged after %d iterations", MAXIT);

278:   if (nrm>1.0) {
279:     alpha = sqrtnrm;
280:     PetscCallCUBLAS(cublasXscal(cublasv2handle,N,&alpha,d_A,one));
281:   }
282:   PetscCall(PetscLogGpuTimeEnd());

284:   PetscCallCUDA(cudaFree(d_M));
285:   PetscCallCUDA(cudaFree(d_M2));
286:   PetscCallCUDA(cudaFree(d_G));
287:   PetscCallCUDA(cudaFree(d_work));
288:   PetscCall(PetscFree(piv));
289:   PetscFunctionReturn(PETSC_SUCCESS);
290: }
291: #endif /* PETSC_HAVE_MAGMA */
292: #endif /* PETSC_HAVE_CUDA */

294: static PetscErrorCode FNEvaluateFunctionMat_Sqrt_Sadeghi(FN fn,Mat A,Mat B)
295: {
296:   PetscBLASInt   n=0;
297:   PetscScalar    *Ba;
298:   PetscInt       m;

300:   PetscFunctionBegin;
301:   if (A!=B) PetscCall(MatCopy(A,B,SAME_NONZERO_PATTERN));
302:   PetscCall(MatDenseGetArray(B,&Ba));
303:   PetscCall(MatGetSize(A,&m,NULL));
304:   PetscCall(PetscBLASIntCast(m,&n));
305:   PetscCall(FNSqrtmSadeghi(fn,n,Ba,n));
306:   PetscCall(MatDenseRestoreArray(B,&Ba));
307:   PetscFunctionReturn(PETSC_SUCCESS);
308: }

310: #if PetscDefined(HAVE_CUDA)
311: PetscErrorCode FNEvaluateFunctionMat_Sqrt_NS_CUDA(FN fn,Mat A,Mat B)
312: {
313:   PetscBLASInt   n=0;
314:   PetscScalar    *Ba;
315:   PetscInt       m;

317:   PetscFunctionBegin;
318:   if (A!=B) PetscCall(MatCopy(A,B,SAME_NONZERO_PATTERN));
319:   PetscCall(MatDenseCUDAGetArray(B,&Ba));
320:   PetscCall(MatGetSize(A,&m,NULL));
321:   PetscCall(PetscBLASIntCast(m,&n));
322:   PetscCall(FNSqrtmNewtonSchulz_CUDA(fn,n,Ba,n,PETSC_FALSE));
323:   PetscCall(MatDenseCUDARestoreArray(B,&Ba));
324:   PetscFunctionReturn(PETSC_SUCCESS);
325: }

327: #if PetscDefined(HAVE_MAGMA)
328: PetscErrorCode FNEvaluateFunctionMat_Sqrt_DBP_CUDAm(FN fn,Mat A,Mat B)
329: {
330:   PetscBLASInt   n=0;
331:   PetscScalar    *T;
332:   PetscInt       m;

334:   PetscFunctionBegin;
335:   if (A!=B) PetscCall(MatCopy(A,B,SAME_NONZERO_PATTERN));
336:   PetscCall(MatDenseCUDAGetArray(B,&T));
337:   PetscCall(MatGetSize(A,&m,NULL));
338:   PetscCall(PetscBLASIntCast(m,&n));
339:   PetscCall(FNSqrtmDenmanBeavers_CUDAm(fn,n,T,n,PETSC_FALSE));
340:   PetscCall(MatDenseCUDARestoreArray(B,&T));
341:   PetscFunctionReturn(PETSC_SUCCESS);
342: }

344: PetscErrorCode FNEvaluateFunctionMat_Sqrt_Sadeghi_CUDAm(FN fn,Mat A,Mat B)
345: {
346:   PetscBLASInt   n=0;
347:   PetscScalar    *Ba;
348:   PetscInt       m;

350:   PetscFunctionBegin;
351:   if (A!=B) PetscCall(MatCopy(A,B,SAME_NONZERO_PATTERN));
352:   PetscCall(MatDenseCUDAGetArray(B,&Ba));
353:   PetscCall(MatGetSize(A,&m,NULL));
354:   PetscCall(PetscBLASIntCast(m,&n));
355:   PetscCall(FNSqrtmSadeghi_CUDAm(fn,n,Ba,n));
356:   PetscCall(MatDenseCUDARestoreArray(B,&Ba));
357:   PetscFunctionReturn(PETSC_SUCCESS);
358: }
359: #endif /* PETSC_HAVE_MAGMA */
360: #endif /* PETSC_HAVE_CUDA */

362: static PetscErrorCode FNView_Sqrt(FN fn,PetscViewer viewer)
363: {
364:   PetscBool      isascii;
365:   char           str[50];
366:   const char     *methodname[] = {
367:                   "Schur method for the square root",
368:                   "Denman-Beavers (product form)",
369:                   "Newton-Schulz iteration",
370:                   "Sadeghi iteration"
371:   };
372:   const int      nmeth=PETSC_STATIC_ARRAY_LENGTH(methodname);

374:   PetscFunctionBegin;
375:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii));
376:   if (isascii) {
377:     if (fn->beta==(PetscScalar)1.0) {
378:       if (fn->alpha==(PetscScalar)1.0) PetscCall(PetscViewerASCIIPrintf(viewer,"  square root: sqrt(x)\n"));
379:       else {
380:         PetscCall(SlepcSNPrintfScalar(str,sizeof(str),fn->alpha,PETSC_TRUE));
381:         PetscCall(PetscViewerASCIIPrintf(viewer,"  square root: sqrt(%s*x)\n",str));
382:       }
383:     } else {
384:       PetscCall(SlepcSNPrintfScalar(str,sizeof(str),fn->beta,PETSC_TRUE));
385:       if (fn->alpha==(PetscScalar)1.0) PetscCall(PetscViewerASCIIPrintf(viewer,"  square root: %s*sqrt(x)\n",str));
386:       else {
387:         PetscCall(PetscViewerASCIIPrintf(viewer,"  square root: %s",str));
388:         PetscCall(PetscViewerASCIIUseTabs(viewer,PETSC_FALSE));
389:         PetscCall(SlepcSNPrintfScalar(str,sizeof(str),fn->alpha,PETSC_TRUE));
390:         PetscCall(PetscViewerASCIIPrintf(viewer,"*sqrt(%s*x)\n",str));
391:         PetscCall(PetscViewerASCIIUseTabs(viewer,PETSC_TRUE));
392:       }
393:     }
394:     if (fn->method<nmeth) PetscCall(PetscViewerASCIIPrintf(viewer,"  computing matrix functions with: %s\n",methodname[fn->method]));
395:   }
396:   PetscFunctionReturn(PETSC_SUCCESS);
397: }

399: /*MC
400:    FNSQRT - FNSQRT = "sqrt" - The square root function $f(x)=\sqrt{x}$.

402:    Level: beginner

404: .seealso: [](sec:fn), `FN`, `FNType`, `FNSetType()`
405: M*/

407: SLEPC_EXTERN PetscErrorCode FNCreate_Sqrt(FN fn)
408: {
409:   PetscFunctionBegin;
410:   fn->ops->evaluatefunction          = FNEvaluateFunction_Sqrt;
411:   fn->ops->evaluatederivative        = FNEvaluateDerivative_Sqrt;
412:   fn->ops->evaluatefunctionmat[0]    = FNEvaluateFunctionMat_Sqrt_Schur;
413:   fn->ops->evaluatefunctionmat[1]    = FNEvaluateFunctionMat_Sqrt_DBP;
414:   fn->ops->evaluatefunctionmat[2]    = FNEvaluateFunctionMat_Sqrt_NS;
415:   fn->ops->evaluatefunctionmat[3]    = FNEvaluateFunctionMat_Sqrt_Sadeghi;
416: #if PetscDefined(HAVE_CUDA)
417:   fn->ops->evaluatefunctionmatcuda[2] = FNEvaluateFunctionMat_Sqrt_NS_CUDA;
418: #if PetscDefined(HAVE_MAGMA)
419:   fn->ops->evaluatefunctionmatcuda[1] = FNEvaluateFunctionMat_Sqrt_DBP_CUDAm;
420:   fn->ops->evaluatefunctionmatcuda[3] = FNEvaluateFunctionMat_Sqrt_Sadeghi_CUDAm;
421: #endif /* PETSC_HAVE_MAGMA */
422: #endif /* PETSC_HAVE_CUDA */
423:   fn->ops->evaluatefunctionmatvec[0] = FNEvaluateFunctionMatVec_Sqrt_Schur;
424:   fn->ops->view                      = FNView_Sqrt;
425:   PetscFunctionReturn(PETSC_SUCCESS);
426: }