| Line | Branch | Exec | Source |
|---|---|---|---|
| 1 | /* | ||
| 2 | - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - | ||
| 3 | SLEPc - Scalable Library for Eigenvalue Problem Computations | ||
| 4 | Copyright (c) 2002-, Universitat Politecnica de Valencia, Spain | ||
| 5 | |||
| 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 | */ | ||
| 13 | |||
| 14 | #include <slepc/private/fnimpl.h> /*I "slepcfn.h" I*/ | ||
| 15 | #include <slepcblaslapack.h> | ||
| 16 | |||
| 17 | 12397 | static PetscErrorCode FNEvaluateFunction_Sqrt(FN fn,PetscScalar x,PetscScalar *y) | |
| 18 | { | ||
| 19 |
1/2✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
|
12397 | PetscFunctionBegin; |
| 20 | #if !defined(PETSC_USE_COMPLEX) | ||
| 21 |
1/4✗ Branch 0 not taken.
✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
|
5033 | PetscCheck(x>=0.0,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Function not defined in the requested value"); |
| 22 | #endif | ||
| 23 | 12397 | *y = PetscSqrtScalar(x); | |
| 24 |
6/12✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 2 times.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 2 times.
✗ Branch 7 not taken.
✓ Branch 8 taken 2 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✓ Branch 11 taken 2 times.
|
12397 | PetscFunctionReturn(PETSC_SUCCESS); |
| 25 | } | ||
| 26 | |||
| 27 | 754 | static PetscErrorCode FNEvaluateDerivative_Sqrt(FN fn,PetscScalar x,PetscScalar *y) | |
| 28 | { | ||
| 29 |
1/2✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
|
754 | PetscFunctionBegin; |
| 30 |
1/4✗ Branch 0 not taken.
✓ Branch 1 taken 10 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
|
754 | PetscCheck(x!=0.0,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Derivative not defined in the requested value"); |
| 31 | #if !defined(PETSC_USE_COMPLEX) | ||
| 32 |
1/4✗ Branch 0 not taken.
✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
|
168 | PetscCheck(x>0.0,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Derivative not defined in the requested value"); |
| 33 | #endif | ||
| 34 | 754 | *y = 1.0/(2.0*PetscSqrtScalar(x)); | |
| 35 |
6/12✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 2 times.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 2 times.
✗ Branch 7 not taken.
✓ Branch 8 taken 2 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✓ Branch 11 taken 2 times.
|
754 | PetscFunctionReturn(PETSC_SUCCESS); |
| 36 | } | ||
| 37 | |||
| 38 | 236 | static PetscErrorCode FNEvaluateFunctionMat_Sqrt_Schur(FN fn,Mat A,Mat B) | |
| 39 | { | ||
| 40 | 236 | PetscBLASInt n=0; | |
| 41 | 236 | PetscScalar *T; | |
| 42 | 236 | PetscInt m; | |
| 43 | |||
| 44 |
1/2✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
|
236 | PetscFunctionBegin; |
| 45 |
5/8✓ Branch 0 taken 10 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 2 times.
✓ Branch 3 taken 8 times.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 2 times.
|
236 | if (A!=B) PetscCall(MatCopy(A,B,SAME_NONZERO_PATTERN)); |
| 46 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
236 | PetscCall(MatDenseGetArray(B,&T)); |
| 47 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
236 | PetscCall(MatGetSize(A,&m,NULL)); |
| 48 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
236 | PetscCall(PetscBLASIntCast(m,&n)); |
| 49 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
236 | PetscCall(FNSqrtmSchur(fn,n,T,n,PETSC_FALSE)); |
| 50 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
236 | PetscCall(MatDenseRestoreArray(B,&T)); |
| 51 |
6/12✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 2 times.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 2 times.
✗ Branch 7 not taken.
✓ Branch 8 taken 2 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✓ Branch 11 taken 2 times.
|
46 | PetscFunctionReturn(PETSC_SUCCESS); |
| 52 | } | ||
| 53 | |||
| 54 | 178 | static PetscErrorCode FNEvaluateFunctionMatVec_Sqrt_Schur(FN fn,Mat A,Vec v) | |
| 55 | { | ||
| 56 | 178 | PetscBLASInt n=0; | |
| 57 | 178 | PetscScalar *T; | |
| 58 | 178 | PetscInt m; | |
| 59 | 178 | Mat B; | |
| 60 | |||
| 61 |
1/2✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
|
178 | PetscFunctionBegin; |
| 62 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
178 | PetscCall(FN_AllocateWorkMat(fn,A,&B)); |
| 63 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
178 | PetscCall(MatDenseGetArray(B,&T)); |
| 64 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
178 | PetscCall(MatGetSize(A,&m,NULL)); |
| 65 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
178 | PetscCall(PetscBLASIntCast(m,&n)); |
| 66 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
178 | PetscCall(FNSqrtmSchur(fn,n,T,n,PETSC_TRUE)); |
| 67 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
178 | PetscCall(MatDenseRestoreArray(B,&T)); |
| 68 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
178 | PetscCall(MatGetColumnVector(B,v,0)); |
| 69 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
178 | PetscCall(FN_FreeWorkMat(fn,&B)); |
| 70 |
6/12✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 2 times.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 2 times.
✗ Branch 7 not taken.
✓ Branch 8 taken 2 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✓ Branch 11 taken 2 times.
|
32 | PetscFunctionReturn(PETSC_SUCCESS); |
| 71 | } | ||
| 72 | |||
| 73 | 132 | static PetscErrorCode FNEvaluateFunctionMat_Sqrt_DBP(FN fn,Mat A,Mat B) | |
| 74 | { | ||
| 75 | 132 | PetscBLASInt n=0; | |
| 76 | 132 | PetscScalar *T; | |
| 77 | 132 | PetscInt m; | |
| 78 | |||
| 79 |
1/2✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
|
132 | PetscFunctionBegin; |
| 80 |
5/8✓ Branch 0 taken 10 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 2 times.
✓ Branch 3 taken 8 times.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 2 times.
|
132 | if (A!=B) PetscCall(MatCopy(A,B,SAME_NONZERO_PATTERN)); |
| 81 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
132 | PetscCall(MatDenseGetArray(B,&T)); |
| 82 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
132 | PetscCall(MatGetSize(A,&m,NULL)); |
| 83 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
132 | PetscCall(PetscBLASIntCast(m,&n)); |
| 84 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
132 | PetscCall(FNSqrtmDenmanBeavers(fn,n,T,n,PETSC_FALSE)); |
| 85 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
132 | PetscCall(MatDenseRestoreArray(B,&T)); |
| 86 |
6/12✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 2 times.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 2 times.
✗ Branch 7 not taken.
✓ Branch 8 taken 2 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✓ Branch 11 taken 2 times.
|
24 | PetscFunctionReturn(PETSC_SUCCESS); |
| 87 | } | ||
| 88 | |||
| 89 | 132 | static PetscErrorCode FNEvaluateFunctionMat_Sqrt_NS(FN fn,Mat A,Mat B) | |
| 90 | { | ||
| 91 | 132 | PetscBLASInt n=0; | |
| 92 | 132 | PetscScalar *Ba; | |
| 93 | 132 | PetscInt m; | |
| 94 | |||
| 95 |
1/2✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
|
132 | PetscFunctionBegin; |
| 96 |
5/8✓ Branch 0 taken 10 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 2 times.
✓ Branch 3 taken 8 times.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 2 times.
|
132 | if (A!=B) PetscCall(MatCopy(A,B,SAME_NONZERO_PATTERN)); |
| 97 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
132 | PetscCall(MatDenseGetArray(B,&Ba)); |
| 98 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
132 | PetscCall(MatGetSize(A,&m,NULL)); |
| 99 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
132 | PetscCall(PetscBLASIntCast(m,&n)); |
| 100 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
132 | PetscCall(FNSqrtmNewtonSchulz(fn,n,Ba,n,PETSC_FALSE)); |
| 101 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
132 | PetscCall(MatDenseRestoreArray(B,&Ba)); |
| 102 |
6/12✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 2 times.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 2 times.
✗ Branch 7 not taken.
✓ Branch 8 taken 2 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✓ Branch 11 taken 2 times.
|
24 | PetscFunctionReturn(PETSC_SUCCESS); |
| 103 | } | ||
| 104 | |||
| 105 | #define MAXIT 50 | ||
| 106 | |||
| 107 | /* | ||
| 108 | Computes the principal square root of the matrix A using the | ||
| 109 | Sadeghi iteration. A is overwritten with sqrtm(A). | ||
| 110 | */ | ||
| 111 | 252 | PetscErrorCode FNSqrtmSadeghi(FN fn,PetscBLASInt n,PetscScalar *A,PetscBLASInt ld) | |
| 112 | { | ||
| 113 | 252 | PetscScalar *M,*M2,*G,*X=A,*work,work1,sqrtnrm; | |
| 114 | 252 | PetscScalar szero=0.0,sone=1.0,smfive=-5.0,s1d16=1.0/16.0; | |
| 115 | 252 | PetscReal tol,Mres=0.0,nrm,rwork[1],done=1.0; | |
| 116 | 252 | PetscInt i,it; | |
| 117 | 252 | PetscBLASInt N,*piv=NULL,lwork=0,query=-1,one=1,zero=0; | |
| 118 | 252 | PetscBool converged=PETSC_FALSE; | |
| 119 | 252 | unsigned int ftz; | |
| 120 | |||
| 121 |
1/2✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
|
252 | PetscFunctionBegin; |
| 122 | 252 | N = n*n; | |
| 123 | 252 | tol = PetscSqrtReal((PetscReal)n)*PETSC_MACHINE_EPSILON/2; | |
| 124 |
3/6✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
252 | PetscCall(SlepcSetFlushToZero(&ftz)); |
| 125 | |||
| 126 | /* query work size */ | ||
| 127 |
13/28✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✓ Branch 3 taken 8 times.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 10 times.
✓ Branch 8 taken 2 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✓ Branch 11 taken 2 times.
✓ Branch 12 taken 2 times.
✗ Branch 13 not taken.
✓ Branch 14 taken 2 times.
✗ Branch 15 not taken.
✗ Branch 16 not taken.
✓ Branch 17 taken 2 times.
✗ Branch 18 not taken.
✗ Branch 19 not taken.
✗ Branch 20 not taken.
✓ Branch 21 taken 2 times.
✗ Branch 22 not taken.
✗ Branch 23 not taken.
✗ Branch 24 not taken.
✓ Branch 25 taken 2 times.
✗ Branch 26 not taken.
✗ Branch 27 not taken.
|
252 | PetscCallLAPACKInfo("LAPACKgetri",LAPACKgetri_(&n,A,&ld,piv,&work1,&query,&info)); |
| 128 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
252 | PetscCall(PetscBLASIntCast((PetscInt)PetscRealPart(work1),&lwork)); |
| 129 | |||
| 130 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
252 | PetscCall(PetscMalloc5(N,&M,N,&M2,N,&G,lwork,&work,n,&piv)); |
| 131 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
252 | PetscCall(PetscArraycpy(M,A,N)); |
| 132 | |||
| 133 | /* scale M */ | ||
| 134 | 252 | nrm = LAPACKlange_("fro",&n,&n,M,&n,rwork); | |
| 135 |
1/2✓ Branch 0 taken 10 times.
✗ Branch 1 not taken.
|
252 | if (nrm>1.0) { |
| 136 | 252 | sqrtnrm = PetscSqrtReal(nrm); | |
| 137 |
13/28✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✓ Branch 3 taken 8 times.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 10 times.
✓ Branch 8 taken 2 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✓ Branch 11 taken 2 times.
✓ Branch 12 taken 2 times.
✗ Branch 13 not taken.
✓ Branch 14 taken 2 times.
✗ Branch 15 not taken.
✗ Branch 16 not taken.
✓ Branch 17 taken 2 times.
✗ Branch 18 not taken.
✗ Branch 19 not taken.
✗ Branch 20 not taken.
✓ Branch 21 taken 2 times.
✗ Branch 22 not taken.
✗ Branch 23 not taken.
✗ Branch 24 not taken.
✓ Branch 25 taken 2 times.
✗ Branch 26 not taken.
✗ Branch 27 not taken.
|
252 | PetscCallLAPACKInfo("LAPACKlascl",LAPACKlascl_("G",&zero,&zero,&nrm,&done,&N,&one,M,&N,&info)); |
| 138 | 252 | tol *= nrm; | |
| 139 | } | ||
| 140 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
252 | PetscCall(PetscInfo(fn,"||A||_F = %g, new tol: %g\n",(double)nrm,(double)tol)); |
| 141 | |||
| 142 | /* X = I */ | ||
| 143 |
4/6✓ Branch 0 taken 5 times.
✓ Branch 1 taken 5 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
252 | PetscCall(PetscArrayzero(X,N)); |
| 144 |
2/2✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
|
13680 | for (i=0;i<n;i++) X[i+i*ld] = 1.0; |
| 145 | |||
| 146 |
2/2✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
|
1328 | for (it=0;it<MAXIT && !converged;it++) { |
| 147 | |||
| 148 | /* G = (5/16)*I + (1/16)*M*(15*I-5*M+M*M) */ | ||
| 149 |
10/20✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 2 times.
✓ Branch 8 taken 2 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✓ Branch 11 taken 2 times.
✓ Branch 12 taken 2 times.
✗ Branch 13 not taken.
✓ Branch 14 taken 2 times.
✗ Branch 15 not taken.
✗ Branch 16 not taken.
✓ Branch 17 taken 2 times.
✗ Branch 18 not taken.
✗ Branch 19 not taken.
|
1076 | PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&n,&n,&n,&sone,M,&ld,M,&ld,&szero,M2,&ld)); |
| 150 |
10/20✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 2 times.
✓ Branch 8 taken 2 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✓ Branch 11 taken 2 times.
✓ Branch 12 taken 2 times.
✗ Branch 13 not taken.
✓ Branch 14 taken 2 times.
✗ Branch 15 not taken.
✗ Branch 16 not taken.
✓ Branch 17 taken 2 times.
✗ Branch 18 not taken.
✗ Branch 19 not taken.
|
1076 | PetscCallBLAS("BLASaxpy",BLASaxpy_(&N,&smfive,M,&one,M2,&one)); |
| 151 |
2/2✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
|
63892 | for (i=0;i<n;i++) M2[i+i*ld] += 15.0; |
| 152 |
10/20✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 2 times.
✓ Branch 8 taken 2 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✓ Branch 11 taken 2 times.
✓ Branch 12 taken 2 times.
✗ Branch 13 not taken.
✓ Branch 14 taken 2 times.
✗ Branch 15 not taken.
✗ Branch 16 not taken.
✓ Branch 17 taken 2 times.
✗ Branch 18 not taken.
✗ Branch 19 not taken.
|
1076 | PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&n,&n,&n,&s1d16,M,&ld,M2,&ld,&szero,G,&ld)); |
| 153 |
2/2✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
|
63892 | for (i=0;i<n;i++) G[i+i*ld] += 5.0/16.0; |
| 154 | |||
| 155 | /* X = X*G */ | ||
| 156 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
1076 | PetscCall(PetscArraycpy(M2,X,N)); |
| 157 |
10/20✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 2 times.
✓ Branch 8 taken 2 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✓ Branch 11 taken 2 times.
✓ Branch 12 taken 2 times.
✗ Branch 13 not taken.
✓ Branch 14 taken 2 times.
✗ Branch 15 not taken.
✗ Branch 16 not taken.
✓ Branch 17 taken 2 times.
✗ Branch 18 not taken.
✗ Branch 19 not taken.
|
1076 | PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&n,&n,&n,&sone,M2,&ld,G,&ld,&szero,X,&ld)); |
| 158 | |||
| 159 | /* M = M*inv(G*G) */ | ||
| 160 |
10/20✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 2 times.
✓ Branch 8 taken 2 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✓ Branch 11 taken 2 times.
✓ Branch 12 taken 2 times.
✗ Branch 13 not taken.
✓ Branch 14 taken 2 times.
✗ Branch 15 not taken.
✗ Branch 16 not taken.
✓ Branch 17 taken 2 times.
✗ Branch 18 not taken.
✗ Branch 19 not taken.
|
1076 | PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&n,&n,&n,&sone,G,&ld,G,&ld,&szero,M2,&ld)); |
| 161 |
13/28✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✓ Branch 3 taken 8 times.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 10 times.
✓ Branch 8 taken 2 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✓ Branch 11 taken 2 times.
✓ Branch 12 taken 2 times.
✗ Branch 13 not taken.
✓ Branch 14 taken 2 times.
✗ Branch 15 not taken.
✗ Branch 16 not taken.
✓ Branch 17 taken 2 times.
✗ Branch 18 not taken.
✗ Branch 19 not taken.
✗ Branch 20 not taken.
✓ Branch 21 taken 2 times.
✗ Branch 22 not taken.
✗ Branch 23 not taken.
✗ Branch 24 not taken.
✓ Branch 25 taken 2 times.
✗ Branch 26 not taken.
✗ Branch 27 not taken.
|
1076 | PetscCallLAPACKInfo("LAPACKgetrf",LAPACKgetrf_(&n,&n,M2,&ld,piv,&info)); |
| 162 |
13/28✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✓ Branch 3 taken 8 times.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 10 times.
✓ Branch 8 taken 2 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✓ Branch 11 taken 2 times.
✓ Branch 12 taken 2 times.
✗ Branch 13 not taken.
✓ Branch 14 taken 2 times.
✗ Branch 15 not taken.
✗ Branch 16 not taken.
✓ Branch 17 taken 2 times.
✗ Branch 18 not taken.
✗ Branch 19 not taken.
✗ Branch 20 not taken.
✓ Branch 21 taken 2 times.
✗ Branch 22 not taken.
✗ Branch 23 not taken.
✗ Branch 24 not taken.
✓ Branch 25 taken 2 times.
✗ Branch 26 not taken.
✗ Branch 27 not taken.
|
1076 | PetscCallLAPACKInfo("LAPACKgetri",LAPACKgetri_(&n,M2,&ld,piv,work,&lwork,&info)); |
| 163 | |||
| 164 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
1076 | PetscCall(PetscArraycpy(G,M,N)); |
| 165 |
10/20✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 2 times.
✓ Branch 8 taken 2 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✓ Branch 11 taken 2 times.
✓ Branch 12 taken 2 times.
✗ Branch 13 not taken.
✓ Branch 14 taken 2 times.
✗ Branch 15 not taken.
✗ Branch 16 not taken.
✓ Branch 17 taken 2 times.
✗ Branch 18 not taken.
✗ Branch 19 not taken.
|
1076 | PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&n,&n,&n,&sone,G,&ld,M2,&ld,&szero,M,&ld)); |
| 166 | |||
| 167 | /* check ||I-M|| */ | ||
| 168 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
1076 | PetscCall(PetscArraycpy(M2,M,N)); |
| 169 |
2/2✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
|
63892 | for (i=0;i<n;i++) M2[i+i*ld] -= 1.0; |
| 170 | 1076 | Mres = LAPACKlange_("fro",&n,&n,M2,&n,rwork); | |
| 171 |
1/4✗ Branch 0 not taken.
✓ Branch 1 taken 10 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
|
1076 | PetscCheck(!PetscIsNanReal(Mres),PETSC_COMM_SELF,PETSC_ERR_FP,"The computed norm is not-a-number"); |
| 172 |
2/2✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
|
1076 | if (Mres<=tol) converged = PETSC_TRUE; |
| 173 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
1076 | PetscCall(PetscInfo(fn,"it: %" PetscInt_FMT " res: %g\n",it,(double)Mres)); |
| 174 |
3/6✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
1076 | 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 | } | ||
| 176 | |||
| 177 |
1/4✗ Branch 0 not taken.
✓ Branch 1 taken 10 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
|
252 | PetscCheck(Mres<=tol,PETSC_COMM_SELF,PETSC_ERR_LIB,"SQRTM not converged after %d iterations",MAXIT); |
| 178 | |||
| 179 | /* undo scaling */ | ||
| 180 |
11/22✓ Branch 0 taken 10 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 2 times.
✓ Branch 3 taken 8 times.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 2 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 2 times.
✓ Branch 10 taken 2 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✓ Branch 13 taken 2 times.
✓ Branch 14 taken 2 times.
✗ Branch 15 not taken.
✓ Branch 16 taken 2 times.
✗ Branch 17 not taken.
✗ Branch 18 not taken.
✓ Branch 19 taken 2 times.
✗ Branch 20 not taken.
✗ Branch 21 not taken.
|
252 | if (nrm>1.0) PetscCallBLAS("BLASscal",BLASscal_(&N,&sqrtnrm,A,&one)); |
| 181 | |||
| 182 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
252 | PetscCall(PetscFree5(M,M2,G,work,piv)); |
| 183 |
3/6✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
74 | PetscCall(SlepcResetFlushToZero(&ftz)); |
| 184 |
6/12✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 2 times.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 2 times.
✗ Branch 7 not taken.
✓ Branch 8 taken 2 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✓ Branch 11 taken 2 times.
|
48 | PetscFunctionReturn(PETSC_SUCCESS); |
| 185 | } | ||
| 186 | |||
| 187 | #if defined(PETSC_HAVE_CUDA) | ||
| 188 | #include "../src/sys/classes/fn/impls/cuda/fnutilcuda.h" | ||
| 189 | #include <slepccupmblas.h> | ||
| 190 | |||
| 191 | #if defined(PETSC_HAVE_MAGMA) | ||
| 192 | #include <slepcmagma.h> | ||
| 193 | |||
| 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 | 48 | PetscErrorCode FNSqrtmSadeghi_CUDAm(FN fn,PetscBLASInt n,PetscScalar *d_A,PetscBLASInt ld) | |
| 200 | { | ||
| 201 | 48 | PetscScalar *d_M,*d_M2,*d_G,*d_work,alpha; | |
| 202 | 48 | const PetscScalar szero=0.0,sone=1.0,smfive=-5.0,s15=15.0,s1d16=1.0/16.0; | |
| 203 | 48 | PetscReal tol,Mres=0.0,nrm,sqrtnrm=1.0; | |
| 204 | 48 | PetscInt it,nb,lwork; | |
| 205 | 48 | PetscBLASInt *piv,N; | |
| 206 | 48 | const PetscBLASInt one=1; | |
| 207 | 48 | PetscBool converged=PETSC_FALSE; | |
| 208 | 48 | cublasHandle_t cublasv2handle; | |
| 209 | |||
| 210 | 48 | PetscFunctionBegin; | |
| 211 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
|
48 | PetscCall(PetscDeviceInitialize(PETSC_DEVICE_CUDA)); /* For CUDA event timers */ |
| 212 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
|
48 | PetscCall(PetscCUBLASGetHandle(&cublasv2handle)); |
| 213 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
|
48 | PetscCall(SlepcMagmaInit()); |
| 214 | 48 | N = n*n; | |
| 215 | 48 | tol = PetscSqrtReal((PetscReal)n)*PETSC_MACHINE_EPSILON/2; | |
| 216 | |||
| 217 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
|
48 | PetscCall(PetscMalloc1(n,&piv)); |
| 218 |
1/4✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
|
48 | PetscCallCUDA(cudaMalloc((void **)&d_M,sizeof(PetscScalar)*N)); |
| 219 |
1/4✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
|
48 | PetscCallCUDA(cudaMalloc((void **)&d_M2,sizeof(PetscScalar)*N)); |
| 220 |
1/4✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
|
48 | PetscCallCUDA(cudaMalloc((void **)&d_G,sizeof(PetscScalar)*N)); |
| 221 | |||
| 222 | 48 | nb = magma_get_xgetri_nb(n); | |
| 223 | 48 | lwork = nb*n; | |
| 224 |
1/4✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
|
48 | PetscCallCUDA(cudaMalloc((void **)&d_work,sizeof(PetscScalar)*lwork)); |
| 225 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
|
48 | PetscCall(PetscLogGpuTimeBegin()); |
| 226 | |||
| 227 | /* M = A */ | ||
| 228 |
1/4✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
|
48 | PetscCallCUDA(cudaMemcpy(d_M,d_A,sizeof(PetscScalar)*N,cudaMemcpyDeviceToDevice)); |
| 229 | |||
| 230 | /* scale M */ | ||
| 231 |
1/10✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
|
48 | PetscCallCUBLAS(cublasXnrm2(cublasv2handle,N,d_M,one,&nrm)); |
| 232 |
1/2✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
|
48 | if (nrm>1.0) { |
| 233 | 48 | sqrtnrm = PetscSqrtReal(nrm); | |
| 234 | 48 | alpha = 1.0/nrm; | |
| 235 |
1/10✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
|
48 | PetscCallCUBLAS(cublasXscal(cublasv2handle,N,&alpha,d_M,one)); |
| 236 | 48 | tol *= nrm; | |
| 237 | } | ||
| 238 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
|
48 | PetscCall(PetscInfo(fn,"||A||_F = %g, new tol: %g\n",(double)nrm,(double)tol)); |
| 239 | |||
| 240 | /* X = I */ | ||
| 241 |
1/4✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
|
48 | PetscCallCUDA(cudaMemset(d_A,0,sizeof(PetscScalar)*N)); |
| 242 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
|
48 | PetscCall(set_diagonal(n,d_A,ld,sone)); |
| 243 | |||
| 244 |
2/2✓ Branch 0 taken 2 times.
✓ Branch 1 taken 2 times.
|
248 | for (it=0;it<MAXIT && !converged;it++) { |
| 245 | |||
| 246 | /* G = (5/16)*I + (1/16)*M*(15*I-5*M+M*M) */ | ||
| 247 |
1/10✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
|
200 | PetscCallCUBLAS(cublasXgemm(cublasv2handle,CUBLAS_OP_N,CUBLAS_OP_N,n,n,n,&sone,d_M,ld,d_M,ld,&szero,d_M2,ld)); |
| 248 |
1/10✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
|
200 | PetscCallCUBLAS(cublasXaxpy(cublasv2handle,N,&smfive,d_M,one,d_M2,one)); |
| 249 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
|
200 | PetscCall(shift_diagonal(n,d_M2,ld,s15)); |
| 250 |
1/10✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
|
200 | PetscCallCUBLAS(cublasXgemm(cublasv2handle,CUBLAS_OP_N,CUBLAS_OP_N,n,n,n,&s1d16,d_M,ld,d_M2,ld,&szero,d_G,ld)); |
| 251 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
|
200 | PetscCall(shift_diagonal(n,d_G,ld,5.0/16.0)); |
| 252 | |||
| 253 | /* X = X*G */ | ||
| 254 |
1/4✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
|
200 | PetscCallCUDA(cudaMemcpy(d_M2,d_A,sizeof(PetscScalar)*N,cudaMemcpyDeviceToDevice)); |
| 255 |
1/10✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
|
200 | PetscCallCUBLAS(cublasXgemm(cublasv2handle,CUBLAS_OP_N,CUBLAS_OP_N,n,n,n,&sone,d_M2,ld,d_G,ld,&szero,d_A,ld)); |
| 256 | |||
| 257 | /* M = M*inv(G*G) */ | ||
| 258 |
1/10✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
|
200 | 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 |
2/6✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 2 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
|
200 | PetscCallMAGMA(magma_xgetrf_gpu,n,n,d_M2,ld,piv); |
| 261 |
2/6✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 2 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
|
200 | PetscCallMAGMA(magma_xgetri_gpu,n,d_M2,ld,piv,d_work,lwork); |
| 262 | /* magma */ | ||
| 263 |
1/4✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
|
200 | PetscCallCUDA(cudaMemcpy(d_G,d_M,sizeof(PetscScalar)*N,cudaMemcpyDeviceToDevice)); |
| 264 |
1/10✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
|
200 | PetscCallCUBLAS(cublasXgemm(cublasv2handle,CUBLAS_OP_N,CUBLAS_OP_N,n,n,n,&sone,d_G,ld,d_M2,ld,&szero,d_M,ld)); |
| 265 | |||
| 266 | /* check ||I-M|| */ | ||
| 267 |
1/4✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
|
200 | PetscCallCUDA(cudaMemcpy(d_M2,d_M,sizeof(PetscScalar)*N,cudaMemcpyDeviceToDevice)); |
| 268 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
|
200 | PetscCall(shift_diagonal(n,d_M2,ld,-1.0)); |
| 269 |
1/10✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
|
200 | PetscCallCUBLAS(cublasXnrm2(cublasv2handle,N,d_M2,one,&Mres)); |
| 270 |
1/4✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
|
200 | PetscCheck(!PetscIsNanReal(Mres),PETSC_COMM_SELF,PETSC_ERR_FP,"The computed norm is not-a-number"); |
| 271 |
2/2✓ Branch 0 taken 2 times.
✓ Branch 1 taken 2 times.
|
200 | if (Mres<=tol) converged = PETSC_TRUE; |
| 272 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
|
200 | PetscCall(PetscInfo(fn,"it: %" PetscInt_FMT " res: %g\n",it,(double)Mres)); |
| 273 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
|
200 | 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 | } | ||
| 275 | |||
| 276 |
1/4✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
|
48 | PetscCheck(Mres<=tol,PETSC_COMM_SELF,PETSC_ERR_LIB,"SQRTM not converged after %d iterations", MAXIT); |
| 277 | |||
| 278 |
1/2✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
|
48 | if (nrm>1.0) { |
| 279 | 48 | alpha = sqrtnrm; | |
| 280 |
1/10✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
|
48 | PetscCallCUBLAS(cublasXscal(cublasv2handle,N,&alpha,d_A,one)); |
| 281 | } | ||
| 282 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
|
48 | PetscCall(PetscLogGpuTimeEnd()); |
| 283 | |||
| 284 |
1/4✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
|
48 | PetscCallCUDA(cudaFree(d_M)); |
| 285 |
1/4✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
|
48 | PetscCallCUDA(cudaFree(d_M2)); |
| 286 |
1/4✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
|
48 | PetscCallCUDA(cudaFree(d_G)); |
| 287 |
1/4✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
|
48 | PetscCallCUDA(cudaFree(d_work)); |
| 288 |
2/4✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 2 times.
|
48 | PetscCall(PetscFree(piv)); |
| 289 | PetscFunctionReturn(PETSC_SUCCESS); | ||
| 290 | } | ||
| 291 | #endif /* PETSC_HAVE_MAGMA */ | ||
| 292 | #endif /* PETSC_HAVE_CUDA */ | ||
| 293 | |||
| 294 | 132 | static PetscErrorCode FNEvaluateFunctionMat_Sqrt_Sadeghi(FN fn,Mat A,Mat B) | |
| 295 | { | ||
| 296 | 132 | PetscBLASInt n=0; | |
| 297 | 132 | PetscScalar *Ba; | |
| 298 | 132 | PetscInt m; | |
| 299 | |||
| 300 |
1/2✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
|
132 | PetscFunctionBegin; |
| 301 |
5/8✓ Branch 0 taken 10 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 2 times.
✓ Branch 3 taken 8 times.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 2 times.
|
132 | if (A!=B) PetscCall(MatCopy(A,B,SAME_NONZERO_PATTERN)); |
| 302 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
132 | PetscCall(MatDenseGetArray(B,&Ba)); |
| 303 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
132 | PetscCall(MatGetSize(A,&m,NULL)); |
| 304 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
132 | PetscCall(PetscBLASIntCast(m,&n)); |
| 305 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
132 | PetscCall(FNSqrtmSadeghi(fn,n,Ba,n)); |
| 306 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
132 | PetscCall(MatDenseRestoreArray(B,&Ba)); |
| 307 |
6/12✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 2 times.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 2 times.
✗ Branch 7 not taken.
✓ Branch 8 taken 2 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✓ Branch 11 taken 2 times.
|
24 | PetscFunctionReturn(PETSC_SUCCESS); |
| 308 | } | ||
| 309 | |||
| 310 | #if defined(PETSC_HAVE_CUDA) | ||
| 311 | 24 | PetscErrorCode FNEvaluateFunctionMat_Sqrt_NS_CUDA(FN fn,Mat A,Mat B) | |
| 312 | { | ||
| 313 | 24 | PetscBLASInt n=0; | |
| 314 | 24 | PetscScalar *Ba; | |
| 315 | 24 | PetscInt m; | |
| 316 | |||
| 317 | 24 | PetscFunctionBegin; | |
| 318 |
2/4✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 2 times.
|
24 | if (A!=B) PetscCall(MatCopy(A,B,SAME_NONZERO_PATTERN)); |
| 319 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
|
24 | PetscCall(MatDenseCUDAGetArray(B,&Ba)); |
| 320 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
|
24 | PetscCall(MatGetSize(A,&m,NULL)); |
| 321 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
|
24 | PetscCall(PetscBLASIntCast(m,&n)); |
| 322 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
|
24 | PetscCall(FNSqrtmNewtonSchulz_CUDA(fn,n,Ba,n,PETSC_FALSE)); |
| 323 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
|
24 | PetscCall(MatDenseCUDARestoreArray(B,&Ba)); |
| 324 | PetscFunctionReturn(PETSC_SUCCESS); | ||
| 325 | } | ||
| 326 | |||
| 327 | #if defined(PETSC_HAVE_MAGMA) | ||
| 328 | 24 | PetscErrorCode FNEvaluateFunctionMat_Sqrt_DBP_CUDAm(FN fn,Mat A,Mat B) | |
| 329 | { | ||
| 330 | 24 | PetscBLASInt n=0; | |
| 331 | 24 | PetscScalar *T; | |
| 332 | 24 | PetscInt m; | |
| 333 | |||
| 334 | 24 | PetscFunctionBegin; | |
| 335 |
2/4✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 2 times.
|
24 | if (A!=B) PetscCall(MatCopy(A,B,SAME_NONZERO_PATTERN)); |
| 336 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
|
24 | PetscCall(MatDenseCUDAGetArray(B,&T)); |
| 337 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
|
24 | PetscCall(MatGetSize(A,&m,NULL)); |
| 338 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
|
24 | PetscCall(PetscBLASIntCast(m,&n)); |
| 339 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
|
24 | PetscCall(FNSqrtmDenmanBeavers_CUDAm(fn,n,T,n,PETSC_FALSE)); |
| 340 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
|
24 | PetscCall(MatDenseCUDARestoreArray(B,&T)); |
| 341 | PetscFunctionReturn(PETSC_SUCCESS); | ||
| 342 | } | ||
| 343 | |||
| 344 | 24 | PetscErrorCode FNEvaluateFunctionMat_Sqrt_Sadeghi_CUDAm(FN fn,Mat A,Mat B) | |
| 345 | { | ||
| 346 | 24 | PetscBLASInt n=0; | |
| 347 | 24 | PetscScalar *Ba; | |
| 348 | 24 | PetscInt m; | |
| 349 | |||
| 350 | 24 | PetscFunctionBegin; | |
| 351 |
2/4✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 2 times.
|
24 | if (A!=B) PetscCall(MatCopy(A,B,SAME_NONZERO_PATTERN)); |
| 352 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
|
24 | PetscCall(MatDenseCUDAGetArray(B,&Ba)); |
| 353 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
|
24 | PetscCall(MatGetSize(A,&m,NULL)); |
| 354 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
|
24 | PetscCall(PetscBLASIntCast(m,&n)); |
| 355 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
|
24 | PetscCall(FNSqrtmSadeghi_CUDAm(fn,n,Ba,n)); |
| 356 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
|
24 | PetscCall(MatDenseCUDARestoreArray(B,&Ba)); |
| 357 | PetscFunctionReturn(PETSC_SUCCESS); | ||
| 358 | } | ||
| 359 | #endif /* PETSC_HAVE_MAGMA */ | ||
| 360 | #endif /* PETSC_HAVE_CUDA */ | ||
| 361 | |||
| 362 | 220 | static PetscErrorCode FNView_Sqrt(FN fn,PetscViewer viewer) | |
| 363 | { | ||
| 364 | 220 | PetscBool isascii; | |
| 365 | 220 | char str[50]; | |
| 366 | 220 | const char *methodname[] = { | |
| 367 | "Schur method for the square root", | ||
| 368 | "Denman-Beavers (product form)", | ||
| 369 | "Newton-Schulz iteration", | ||
| 370 | "Sadeghi iteration" | ||
| 371 | }; | ||
| 372 | 220 | const int nmeth=PETSC_STATIC_ARRAY_LENGTH(methodname); | |
| 373 | |||
| 374 |
1/2✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
|
220 | PetscFunctionBegin; |
| 375 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
220 | PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii)); |
| 376 |
1/2✓ Branch 0 taken 10 times.
✗ Branch 1 not taken.
|
220 | if (isascii) { |
| 377 |
2/2✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
|
220 | if (fn->beta==(PetscScalar)1.0) { |
| 378 |
5/8✓ Branch 0 taken 10 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 2 times.
✓ Branch 3 taken 8 times.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 2 times.
|
20 | 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 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
200 | PetscCall(SlepcSNPrintfScalar(str,sizeof(str),fn->beta,PETSC_TRUE)); |
| 385 |
1/8✗ Branch 0 not taken.
✓ Branch 1 taken 10 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
|
200 | if (fn->alpha==(PetscScalar)1.0) PetscCall(PetscViewerASCIIPrintf(viewer," square root: %s*sqrt(x)\n",str)); |
| 386 | else { | ||
| 387 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
200 | PetscCall(PetscViewerASCIIPrintf(viewer," square root: %s",str)); |
| 388 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
200 | PetscCall(PetscViewerASCIIUseTabs(viewer,PETSC_FALSE)); |
| 389 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
200 | PetscCall(SlepcSNPrintfScalar(str,sizeof(str),fn->alpha,PETSC_TRUE)); |
| 390 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
200 | PetscCall(PetscViewerASCIIPrintf(viewer,"*sqrt(%s*x)\n",str)); |
| 391 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
200 | PetscCall(PetscViewerASCIIUseTabs(viewer,PETSC_TRUE)); |
| 392 | } | ||
| 393 | } | ||
| 394 |
5/8✓ Branch 0 taken 10 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 2 times.
✓ Branch 3 taken 8 times.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 2 times.
|
220 | if (fn->method<nmeth) PetscCall(PetscViewerASCIIPrintf(viewer," computing matrix functions with: %s\n",methodname[fn->method])); |
| 395 | } | ||
| 396 |
6/12✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 2 times.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 2 times.
✗ Branch 7 not taken.
✓ Branch 8 taken 2 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✓ Branch 11 taken 2 times.
|
32 | PetscFunctionReturn(PETSC_SUCCESS); |
| 397 | } | ||
| 398 | |||
| 399 | /*MC | ||
| 400 | FNSQRT - FNSQRT = "sqrt" - The square root function $f(x)=\sqrt{x}$. | ||
| 401 | |||
| 402 | Level: beginner | ||
| 403 | |||
| 404 | .seealso: [](sec:fn), `FN`, `FNType`, `FNSetType()` | ||
| 405 | M*/ | ||
| 406 | |||
| 407 | 404 | SLEPC_EXTERN PetscErrorCode FNCreate_Sqrt(FN fn) | |
| 408 | { | ||
| 409 |
1/2✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
|
404 | PetscFunctionBegin; |
| 410 | 404 | fn->ops->evaluatefunction = FNEvaluateFunction_Sqrt; | |
| 411 | 404 | fn->ops->evaluatederivative = FNEvaluateDerivative_Sqrt; | |
| 412 | 404 | fn->ops->evaluatefunctionmat[0] = FNEvaluateFunctionMat_Sqrt_Schur; | |
| 413 | 404 | fn->ops->evaluatefunctionmat[1] = FNEvaluateFunctionMat_Sqrt_DBP; | |
| 414 | 404 | fn->ops->evaluatefunctionmat[2] = FNEvaluateFunctionMat_Sqrt_NS; | |
| 415 | 404 | fn->ops->evaluatefunctionmat[3] = FNEvaluateFunctionMat_Sqrt_Sadeghi; | |
| 416 | #if defined(PETSC_HAVE_CUDA) | ||
| 417 | 88 | fn->ops->evaluatefunctionmatcuda[2] = FNEvaluateFunctionMat_Sqrt_NS_CUDA; | |
| 418 | #if defined(PETSC_HAVE_MAGMA) | ||
| 419 | 88 | fn->ops->evaluatefunctionmatcuda[1] = FNEvaluateFunctionMat_Sqrt_DBP_CUDAm; | |
| 420 | 88 | fn->ops->evaluatefunctionmatcuda[3] = FNEvaluateFunctionMat_Sqrt_Sadeghi_CUDAm; | |
| 421 | #endif /* PETSC_HAVE_MAGMA */ | ||
| 422 | #endif /* PETSC_HAVE_CUDA */ | ||
| 423 | 404 | fn->ops->evaluatefunctionmatvec[0] = FNEvaluateFunctionMatVec_Sqrt_Schur; | |
| 424 | 404 | fn->ops->view = FNView_Sqrt; | |
| 425 |
6/12✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 2 times.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 2 times.
✗ Branch 7 not taken.
✓ Branch 8 taken 2 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✓ Branch 11 taken 2 times.
|
404 | PetscFunctionReturn(PETSC_SUCCESS); |
| 426 | } | ||
| 427 |