Actual source code: densecuda.cu

petsc-3.12.4 2020-02-04
Report Typos and Errors
  1: /*
  2:      Defines the matrix operations for sequential dense with CUDA
  3: */
  4: #include <petscpkg_version.h>
  5: #define PETSC_SKIP_IMMINTRIN_H_CUDAWORKAROUND 1
  6:  #include <../src/mat/impls/dense/seq/dense.h>
  7:  #include <petsccublas.h>

  9: /* cublas definitions are here */
 10:  #include <../src/vec/vec/impls/seq/seqcuda/cudavecimpl.h>

 12: #if defined(PETSC_USE_COMPLEX)
 13: #if defined(PETSC_USE_REAL_SINGLE)
 14: #define cusolverDnXpotrf(a,b,c,d,e,f,g,h)        cusolverDnCpotrf((a),(b),(c),(cuComplex*)(d),(e),(cuComplex*)(f),(g),(h))
 15: #define cusolverDnXpotrf_bufferSize(a,b,c,d,e,f) cusolverDnCpotrf_bufferSize((a),(b),(c),(cuComplex*)(d),(e),(f))
 16: #define cusolverDnXpotrs(a,b,c,d,e,f,g,h,i)      cusolverDnCpotrs((a),(b),(c),(d),(cuComplex*)(e),(f),(cuComplex*)(g),(h),(i))
 17: #define cusolverDnXpotri(a,b,c,d,e,f,g,h)        cusolverDnCpotri((a),(b),(c),(cuComplex*)(d),(e),(cuComplex*)(f),(g),(h))
 18: #define cusolverDnXpotri_bufferSize(a,b,c,d,e,f) cusolverDnCpotri_bufferSize((a),(b),(c),(cuComplex*)(d),(e),(f))
 19: #define cusolverDnXsytrf(a,b,c,d,e,f,g,h,i)      cusolverDnCsytrf((a),(b),(c),(cuComplex*)(d),(e),(f),(cuComplex*)(g),(h),(i))
 20: #define cusolverDnXsytrf_bufferSize(a,b,c,d,e)   cusolverDnCsytrf_bufferSize((a),(b),(cuComplex*)(c),(d),(e))
 21: #define cusolverDnXgetrf(a,b,c,d,e,f,g,h)        cusolverDnCgetrf((a),(b),(c),(cuComplex*)(d),(e),(cuComplex*)(f),(g),(h))
 22: #define cusolverDnXgetrf_bufferSize(a,b,c,d,e,f) cusolverDnCgetrf_bufferSize((a),(b),(c),(cuComplex*)(d),(e),(f))
 23: #define cusolverDnXgetrs(a,b,c,d,e,f,g,h,i,j)    cusolverDnCgetrs((a),(b),(c),(d),(cuComplex*)(e),(f),(g),(cuComplex*)(h),(i),(j))
 24: #else /* complex double */
 25: #define cusolverDnXpotrf(a,b,c,d,e,f,g,h)        cusolverDnZpotrf((a),(b),(c),(cuDoubleComplex*)(d),(e),(cuDoubleComplex*)(f),(g),(h))
 26: #define cusolverDnXpotrf_bufferSize(a,b,c,d,e,f) cusolverDnZpotrf_bufferSize((a),(b),(c),(cuDoubleComplex*)(d),(e),(f))
 27: #define cusolverDnXpotrs(a,b,c,d,e,f,g,h,i)      cusolverDnZpotrs((a),(b),(c),(d),(cuDoubleComplex*)(e),(f),(cuDoubleComplex*)(g),(h),(i))
 28: #define cusolverDnXpotri(a,b,c,d,e,f,g,h)        cusolverDnZpotri((a),(b),(c),(cuDoubleComplex*)(d),(e),(cuDoubleComplex*)(f),(g),(h))
 29: #define cusolverDnXpotri_bufferSize(a,b,c,d,e,f) cusolverDnZpotri_bufferSize((a),(b),(c),(cuDoubleComplex*)(d),(e),(f))
 30: #define cusolverDnXsytrf(a,b,c,d,e,f,g,h,i)      cusolverDnZsytrf((a),(b),(c),(cuDoubleComplex*)(d),(e),(f),(cuDoubleComplex*)(g),(h),(i))
 31: #define cusolverDnXsytrf_bufferSize(a,b,c,d,e)   cusolverDnZsytrf_bufferSize((a),(b),(cuDoubleComplex*)(c),(d),(e))
 32: #define cusolverDnXgetrf(a,b,c,d,e,f,g,h)        cusolverDnZgetrf((a),(b),(c),(cuDoubleComplex*)(d),(e),(cuDoubleComplex*)(f),(g),(h))
 33: #define cusolverDnXgetrf_bufferSize(a,b,c,d,e,f) cusolverDnZgetrf_bufferSize((a),(b),(c),(cuDoubleComplex*)(d),(e),(f))
 34: #define cusolverDnXgetrs(a,b,c,d,e,f,g,h,i,j)    cusolverDnZgetrs((a),(b),(c),(d),(cuDoubleComplex*)(e),(f),(g),(cuDoubleComplex*)(h),(i),(j))
 35: #endif
 36: #else /* real single */
 37: #if defined(PETSC_USE_REAL_SINGLE)
 38: #define cusolverDnXpotrf(a,b,c,d,e,f,g,h)        cusolverDnSpotrf((a),(b),(c),(d),(e),(f),(g),(h))
 39: #define cusolverDnXpotrf_bufferSize(a,b,c,d,e,f) cusolverDnSpotrf_bufferSize((a),(b),(c),(d),(e),(f))
 40: #define cusolverDnXpotrs(a,b,c,d,e,f,g,h,i)      cusolverDnSpotrs((a),(b),(c),(d),(e),(f),(g),(h),(i))
 41: #define cusolverDnXpotri(a,b,c,d,e,f,g,h)        cusolverDnSpotri((a),(b),(c),(d),(e),(f),(g),(h))
 42: #define cusolverDnXpotri_bufferSize(a,b,c,d,e,f) cusolverDnSpotri_bufferSize((a),(b),(c),(d),(e),(f))
 43: #define cusolverDnXsytrf(a,b,c,d,e,f,g,h,i)      cusolverDnSsytrf((a),(b),(c),(d),(e),(f),(g),(h),(i))
 44: #define cusolverDnXsytrf_bufferSize(a,b,c,d,e)   cusolverDnSsytrf_bufferSize((a),(b),(c),(d),(e))
 45: #define cusolverDnXgetrf(a,b,c,d,e,f,g,h)        cusolverDnSgetrf((a),(b),(c),(d),(e),(f),(g),(h))
 46: #define cusolverDnXgetrf_bufferSize(a,b,c,d,e,f) cusolverDnSgetrf_bufferSize((a),(b),(c),(d),(e),(f))
 47: #define cusolverDnXgetrs(a,b,c,d,e,f,g,h,i,j)    cusolverDnSgetrs((a),(b),(c),(d),(e),(f),(g),(h),(i),(j))
 48: #else /* real double */
 49: #define cusolverDnXpotrf(a,b,c,d,e,f,g,h)        cusolverDnDpotrf((a),(b),(c),(d),(e),(f),(g),(h))
 50: #define cusolverDnXpotrf_bufferSize(a,b,c,d,e,f) cusolverDnDpotrf_bufferSize((a),(b),(c),(d),(e),(f))
 51: #define cusolverDnXpotrs(a,b,c,d,e,f,g,h,i)      cusolverDnDpotrs((a),(b),(c),(d),(e),(f),(g),(h),(i))
 52: #define cusolverDnXpotri(a,b,c,d,e,f,g,h)        cusolverDnDpotri((a),(b),(c),(d),(e),(f),(g),(h))
 53: #define cusolverDnXpotri_bufferSize(a,b,c,d,e,f) cusolverDnDpotri_bufferSize((a),(b),(c),(d),(e),(f))
 54: #define cusolverDnXsytrf(a,b,c,d,e,f,g,h,i)      cusolverDnDsytrf((a),(b),(c),(d),(e),(f),(g),(h),(i))
 55: #define cusolverDnXsytrf_bufferSize(a,b,c,d,e)   cusolverDnDsytrf_bufferSize((a),(b),(c),(d),(e))
 56: #define cusolverDnXgetrf(a,b,c,d,e,f,g,h)        cusolverDnDgetrf((a),(b),(c),(d),(e),(f),(g),(h))
 57: #define cusolverDnXgetrf_bufferSize(a,b,c,d,e,f) cusolverDnDgetrf_bufferSize((a),(b),(c),(d),(e),(f))
 58: #define cusolverDnXgetrs(a,b,c,d,e,f,g,h,i,j)    cusolverDnDgetrs((a),(b),(c),(d),(e),(f),(g),(h),(i),(j))
 59: #endif
 60: #endif

 62: typedef struct {
 63:   PetscScalar *d_v;   /* pointer to the matrix on the GPU */
 64:   /* factorization support */
 65:   int         *d_fact_ipiv; /* device pivots */
 66:   PetscScalar *d_fact_work; /* device workspace */
 67:   int         fact_lwork;
 68:   int         *d_fact_info; /* device info */
 69:   /* workspace */
 70:   Vec         workvec;
 71: } Mat_SeqDenseCUDA;

 73: PetscErrorCode MatSeqDenseCUDACopyFromGPU(Mat A)
 74: {
 75:   Mat_SeqDense     *cA = (Mat_SeqDense*)A->data;
 76:   Mat_SeqDenseCUDA *dA = (Mat_SeqDenseCUDA*)A->spptr;
 77:   PetscErrorCode   ierr;
 78:   cudaError_t      cerr;

 82:   PetscInfo3(A,"%s matrix %d x %d\n",A->offloadmask == PETSC_OFFLOAD_GPU ? "Copy" : "Reusing",A->rmap->n,A->cmap->n);
 83:   if (A->offloadmask == PETSC_OFFLOAD_GPU) {
 84:     PetscLogEventBegin(MAT_DenseCopyFromGPU,A,0,0,0);
 85:     if (cA->lda > A->rmap->n) {
 86:       PetscInt j,m = A->rmap->n;

 88:       for (j=0; j<A->cmap->n; j++) { /* TODO: it can be done better */
 89:         cerr = cudaMemcpy(cA->v + j*cA->lda,dA->d_v + j*cA->lda,m*sizeof(PetscScalar),cudaMemcpyDeviceToHost);CHKERRCUDA(cerr);
 90:       }
 91:     } else {
 92:       cerr = cudaMemcpy(cA->v,dA->d_v,cA->lda*sizeof(PetscScalar)*A->cmap->n,cudaMemcpyDeviceToHost);CHKERRCUDA(cerr);
 93:     }
 94:     PetscLogGpuToCpu(cA->lda*sizeof(PetscScalar)*A->cmap->n);
 95:     PetscLogEventEnd(MAT_DenseCopyFromGPU,A,0,0,0);

 97:     A->offloadmask = PETSC_OFFLOAD_BOTH;
 98:   }
 99:   return(0);
100: }

102: PetscErrorCode MatSeqDenseCUDACopyToGPU(Mat A)
103: {
104:   Mat_SeqDense     *cA = (Mat_SeqDense*)A->data;
105:   Mat_SeqDenseCUDA *dA = (Mat_SeqDenseCUDA*)A->spptr;
106:   PetscBool        copy;
107:   PetscErrorCode   ierr;
108:   cudaError_t      cerr;

112:   if (A->pinnedtocpu) return(0);
113:   if (!dA->d_v) {
114:     cerr = cudaMalloc((void**)&dA->d_v,cA->lda*cA->Nmax*sizeof(PetscScalar));CHKERRCUDA(cerr);
115:   }
116:   copy = (PetscBool)(A->offloadmask == PETSC_OFFLOAD_CPU || A->offloadmask == PETSC_OFFLOAD_UNALLOCATED);
117:   PetscInfo3(A,"%s matrix %d x %d\n",copy ? "Copy" : "Reusing",A->rmap->n,A->cmap->n);
118:   if (copy) {
119:     PetscLogEventBegin(MAT_DenseCopyToGPU,A,0,0,0);
120:     if (cA->lda > A->rmap->n) {
121:       PetscInt j,m = A->rmap->n;

123:       for (j=0; j<A->cmap->n; j++) { /* TODO: it can be done better */
124:         cerr = cudaMemcpy(dA->d_v + j*cA->lda,cA->v + j*cA->lda,m*sizeof(PetscScalar),cudaMemcpyHostToDevice);CHKERRCUDA(cerr);
125:       }
126:     } else {
127:       cerr = cudaMemcpy(dA->d_v,cA->v,cA->lda*sizeof(PetscScalar)*A->cmap->n,cudaMemcpyHostToDevice);CHKERRCUDA(cerr);
128:     }
129:     PetscLogCpuToGpu(cA->lda*sizeof(PetscScalar)*A->cmap->n);
130:     PetscLogEventEnd(MAT_DenseCopyToGPU,A,0,0,0);

132:     A->offloadmask = PETSC_OFFLOAD_BOTH;
133:   }
134:   return(0);
135: }

137: PetscErrorCode MatDenseCUDAGetArrayWrite(Mat A, PetscScalar **a)
138: {
139:   Mat_SeqDenseCUDA *dA = (Mat_SeqDenseCUDA*)A->spptr;

143:   if (!dA->d_v) {
144:     Mat_SeqDense *cA = (Mat_SeqDense*)A->data;
145:     cudaError_t  cerr;

147:     cerr = cudaMalloc((void**)&dA->d_v,cA->lda*cA->Nmax*sizeof(PetscScalar));CHKERRCUDA(cerr);
148:   }
149:   *a = dA->d_v;
150:   return(0);
151: }

153: PetscErrorCode MatDenseCUDARestoreArrayWrite(Mat A, PetscScalar **a)
154: {
156:   *a = NULL;
157:   A->offloadmask = PETSC_OFFLOAD_GPU;
158:   return(0);
159: }

161: PetscErrorCode MatDenseCUDAGetArrayRead(Mat A, const PetscScalar **a)
162: {
163:   Mat_SeqDenseCUDA *dA = (Mat_SeqDenseCUDA*)A->spptr;
164:   PetscErrorCode   ierr;

168:   MatSeqDenseCUDACopyToGPU(A);
169:   *a   = dA->d_v;
170:   return(0);
171: }

173: PetscErrorCode MatDenseCUDARestoreArrayRead(Mat A, const PetscScalar **a)
174: {
176:   *a = NULL;
177:   return(0);
178: }

180: PetscErrorCode MatDenseCUDAGetArray(Mat A, PetscScalar **a)
181: {
182:   Mat_SeqDenseCUDA *dA = (Mat_SeqDenseCUDA*)A->spptr;
183:   PetscErrorCode   ierr;

187:   MatSeqDenseCUDACopyToGPU(A);
188:   *a   = dA->d_v;
189:   return(0);
190: }

192: PetscErrorCode MatDenseCUDARestoreArray(Mat A, PetscScalar **a)
193: {
195:   *a = NULL;
196:   A->offloadmask = PETSC_OFFLOAD_GPU;
197:   return(0);
198: }

200: PETSC_EXTERN PetscErrorCode MatSeqDenseCUDAInvertFactors_Private(Mat A)
201: {
202: #if PETSC_PKG_CUDA_VERSION_GE(10,1,0)
203:   Mat_SeqDense       *a = (Mat_SeqDense*)A->data;
204:   Mat_SeqDenseCUDA   *dA = (Mat_SeqDenseCUDA*)A->spptr;
205:   PetscScalar        *da;
206:   PetscErrorCode     ierr;
207:   cudaError_t        ccer;
208:   cusolverStatus_t   cerr;
209:   cusolverDnHandle_t handle;
210:   int                n,lda;
211: #if defined(PETSC_USE_DEBUG)
212:   int                info;
213: #endif

216:   if (!A->rmap->n || !A->cmap->n) return(0);
217:   PetscCUSOLVERDnGetHandle(&handle);
218:   PetscMPIIntCast(A->cmap->n,&n);
219:   PetscMPIIntCast(a->lda,&lda);
220:   if (A->factortype == MAT_FACTOR_LU) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"cusolverDngetri not implemented");
221:   else if (A->factortype == MAT_FACTOR_CHOLESKY) {
222:     if (!dA->d_fact_ipiv) { /* spd */
223:       int il;

225:       MatDenseCUDAGetArray(A,&da);
226:       cerr = cusolverDnXpotri_bufferSize(handle,CUBLAS_FILL_MODE_LOWER,n,da,lda,&il);CHKERRCUSOLVER(cerr);
227:       if (il > dA->fact_lwork) {
228:         dA->fact_lwork = il;

230:         ccer = cudaFree(dA->d_fact_work);CHKERRCUDA(ccer);
231:         ccer = cudaMalloc((void**)&dA->d_fact_work,dA->fact_lwork*sizeof(*dA->d_fact_work));CHKERRCUDA(ccer);
232:       }
233:       PetscLogGpuTimeBegin();
234:       cerr = cusolverDnXpotri(handle,CUBLAS_FILL_MODE_LOWER,n,da,lda,dA->d_fact_work,dA->fact_lwork,dA->d_fact_info);CHKERRCUSOLVER(cerr);
235:       WaitForGPU();CHKERRCUDA(ierr);
236:       PetscLogGpuTimeEnd();
237:       MatDenseCUDARestoreArray(A,&da);
238:       /* TODO (write cuda kernel) */
239:       MatSeqDenseSymmetrize_Private(A,PETSC_TRUE);
240:     } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"cusolverDnsytri not implemented");
241:   }
242: #if defined(PETSC_USE_DEBUG)
243:   ccer = cudaMemcpy(&info, dA->d_fact_info, sizeof(int), cudaMemcpyDeviceToHost);CHKERRCUDA(ccer);
244:   if (info > 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_MAT_CH_ZRPVT,"Bad factorization: leading minor of order %d is zero",info);
245:   else if (info < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Wrong argument to cuSolver %d",-info);
246: #endif
247:   PetscLogGpuFlops(1.0*n*n*n/3.0);
248:   A->ops->solve          = NULL;
249:   A->ops->solvetranspose = NULL;
250:   A->ops->matsolve       = NULL;
251:   A->factortype          = MAT_FACTOR_NONE;

253:   PetscFree(A->solvertype);
254:   return(0);
255: #else
256:   SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Upgrade to CUDA version 10.1.0 or higher");
257: #endif
258: }

260: static PetscErrorCode MatMatSolve_SeqDenseCUDA(Mat A,Mat B,Mat X)
261: {
262:   Mat_SeqDense       *a = (Mat_SeqDense*)A->data;
263:   Mat_SeqDense       *x = (Mat_SeqDense*)X->data;
264:   Mat_SeqDenseCUDA   *dA = (Mat_SeqDenseCUDA*)A->spptr;
265:   const PetscScalar  *da;
266:   PetscScalar        *dx;
267:   cusolverDnHandle_t handle;
268:   PetscBool          iscuda;
269:   int                nrhs,n,lda,ldx;
270: #if defined(PETSC_USE_DEBUG)
271:   int                info;
272:   cudaError_t        ccer;
273: #endif
274:   cusolverStatus_t   cerr;
275:   PetscErrorCode     ierr;

278:   if (A->factortype == MAT_FACTOR_NONE) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix must be factored to solve");
279:   if (!dA->d_fact_work) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix must be factored to solve");
280:   PetscObjectTypeCompareAny((PetscObject)X,&iscuda,VECSEQCUDA,VECMPICUDA,"");
281:   MatCopy(B,X,SAME_NONZERO_PATTERN);
282:   MatDenseCUDAGetArrayRead(A,&da);
283:   /* MatMatSolve does not have a dispatching mechanism, we may end up with a MATSEQDENSE here */
284:   PetscObjectTypeCompare((PetscObject)X,MATSEQDENSECUDA,&iscuda);
285:   if (!iscuda) {
286:     MatConvert(X,MATSEQDENSECUDA,MAT_INPLACE_MATRIX,&X);
287:   }
288:   MatDenseCUDAGetArray(X,&dx);
289:   PetscMPIIntCast(A->rmap->n,&n);
290:   PetscMPIIntCast(X->cmap->n,&nrhs);
291:   PetscMPIIntCast(a->lda,&lda);
292:   PetscMPIIntCast(x->lda,&ldx);
293:   PetscCUSOLVERDnGetHandle(&handle);
294:   PetscLogGpuTimeBegin();
295:   if (A->factortype == MAT_FACTOR_LU) {
296:     PetscInfo2(A,"LU solve %d x %d on backend\n",n,n);
297:     cerr = cusolverDnXgetrs(handle,CUBLAS_OP_N,n,nrhs,da,lda,dA->d_fact_ipiv,dx,ldx,dA->d_fact_info);CHKERRCUSOLVER(cerr);
298:   } else if (A->factortype == MAT_FACTOR_CHOLESKY) {
299:     PetscInfo2(A,"Cholesky solve %d x %d on backend\n",n,n);
300:     if (!dA->d_fact_ipiv) { /* spd */
301:       /* ========= Program hit cudaErrorNotReady (error 34) due to "device not ready" on CUDA API call to cudaEventQuery. */
302:       cerr = cusolverDnXpotrs(handle,CUBLAS_FILL_MODE_LOWER,n,nrhs,da,lda,dx,ldx,dA->d_fact_info);CHKERRCUSOLVER(cerr);
303:     } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"cusolverDnsytrs not implemented");
304:   } else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Unknown factor type %d",A->factortype);
305:   WaitForGPU();CHKERRCUDA(ierr);
306:   PetscLogGpuTimeEnd();
307:   MatDenseCUDARestoreArrayRead(A,&da);
308:   MatDenseCUDARestoreArray(X,&dx);
309:   if (!iscuda) {
310:     MatConvert(X,MATSEQDENSE,MAT_INPLACE_MATRIX,&X);
311:   }
312: #if defined(PETSC_USE_DEBUG)
313:   ccer = cudaMemcpy(&info, dA->d_fact_info, sizeof(int), cudaMemcpyDeviceToHost);CHKERRCUDA(ccer);
314:   if (info > 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_MAT_CH_ZRPVT,"Bad factorization: zero pivot in row %d",info-1);
315:   else if (info < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Wrong argument to cuSolver %d",-info);
316: #endif
317:   PetscLogGpuFlops(nrhs*(2.0*n*n - n));
318:   return(0);
319: }

321: static PetscErrorCode MatSolve_SeqDenseCUDA_Private(Mat A,Vec xx,Vec yy,PetscBool trans)
322: {
323:   Mat_SeqDense       *a = (Mat_SeqDense*)A->data;
324:   Mat_SeqDenseCUDA   *dA = (Mat_SeqDenseCUDA*)A->spptr;
325:   const PetscScalar  *da;
326:   PetscScalar        *y;
327:   cusolverDnHandle_t handle;
328:   int                one = 1,n,lda;
329: #if defined(PETSC_USE_DEBUG)
330:   int                info;
331:   cudaError_t        ccer;
332: #endif
333:   cusolverStatus_t   cerr;
334:   PetscBool          iscuda;
335:   PetscErrorCode     ierr;

338:   if (A->factortype == MAT_FACTOR_NONE) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix must be factored to solve");
339:   if (!dA->d_fact_work) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix must be factored to solve");
340:   PetscMPIIntCast(A->rmap->n,&n);
341:   /* MatSolve does not have a dispatching mechanism, we may end up with a VECSTANDARD here */
342:   PetscObjectTypeCompareAny((PetscObject)yy,&iscuda,VECSEQCUDA,VECMPICUDA,"");
343:   if (iscuda) {
344:     VecCopy(xx,yy);
345:     VecCUDAGetArray(yy,&y);
346:   } else {
347:     if (!dA->workvec) {
348:       MatCreateVecs(A,&dA->workvec,NULL);
349:     }
350:     VecCopy(xx,dA->workvec);
351:     VecCUDAGetArray(dA->workvec,&y);
352:   }
353:   MatDenseCUDAGetArrayRead(A,&da);
354:   PetscMPIIntCast(a->lda,&lda);
355:   PetscCUSOLVERDnGetHandle(&handle);
356:   PetscLogGpuTimeBegin();
357:   if (A->factortype == MAT_FACTOR_LU) {
358:     PetscInfo2(A,"LU solve %d x %d on backend\n",n,n);
359:     cerr = cusolverDnXgetrs(handle,trans ? CUBLAS_OP_T : CUBLAS_OP_N,n,one,da,lda,dA->d_fact_ipiv,y,n,dA->d_fact_info);CHKERRCUSOLVER(cerr);
360:   } else if (A->factortype == MAT_FACTOR_CHOLESKY) {
361:     PetscInfo2(A,"Cholesky solve %d x %d on backend\n",n,n);
362:     if (!dA->d_fact_ipiv) { /* spd */
363:       /* ========= Program hit cudaErrorNotReady (error 34) due to "device not ready" on CUDA API call to cudaEventQuery. */
364:       cerr = cusolverDnXpotrs(handle,CUBLAS_FILL_MODE_LOWER,n,one,da,lda,y,n,dA->d_fact_info);CHKERRCUSOLVER(cerr);
365:     } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"cusolverDnsytrs not implemented");
366:   } else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Unknown factor type %d",A->factortype);
367:   WaitForGPU();CHKERRCUDA(ierr);
368:   PetscLogGpuTimeEnd();
369:   if (iscuda) {
370:     VecCUDARestoreArray(yy,&y);
371:   } else {
372:     VecCUDARestoreArray(dA->workvec,&y);
373:     VecCopy(dA->workvec,yy);
374:   }
375:   MatDenseCUDARestoreArrayRead(A,&da);
376: #if defined(PETSC_USE_DEBUG)
377:   ccer = cudaMemcpy(&info, dA->d_fact_info, sizeof(int), cudaMemcpyDeviceToHost);CHKERRCUDA(ccer);
378:   if (info > 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_MAT_CH_ZRPVT,"Bad factorization: zero pivot in row %d",info-1);
379:   else if (info < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Wrong argument to cuSolver %d",-info);
380: #endif
381:   PetscLogGpuFlops(2.0*n*n - n);
382:   return(0);
383: }

385: static PetscErrorCode MatSolve_SeqDenseCUDA(Mat A,Vec xx,Vec yy)
386: {
387:   PetscErrorCode     ierr;

390:   MatSolve_SeqDenseCUDA_Private(A,xx,yy,PETSC_FALSE);
391:   return(0);
392: }

394: static PetscErrorCode MatSolveTranspose_SeqDenseCUDA(Mat A,Vec xx,Vec yy)
395: {
396:   PetscErrorCode     ierr;

399:   MatSolve_SeqDenseCUDA_Private(A,xx,yy,PETSC_TRUE);
400:   return(0);
401: }

403: static PetscErrorCode MatLUFactor_SeqDenseCUDA(Mat A,IS rperm,IS cperm,const MatFactorInfo *factinfo)
404: {
405:   Mat_SeqDense       *a = (Mat_SeqDense*)A->data;
406:   Mat_SeqDenseCUDA   *dA = (Mat_SeqDenseCUDA*)A->spptr;
407:   PetscScalar        *da;
408:   int                m,n,lda;
409: #if defined(PETSC_USE_DEBUG)
410:   int                info;
411: #endif
412:   cusolverStatus_t   cerr;
413:   cusolverDnHandle_t handle;
414:   cudaError_t        ccer;
415:   PetscErrorCode     ierr;

418:   if (!A->rmap->n || !A->cmap->n) return(0);
419:   PetscCUSOLVERDnGetHandle(&handle);
420:   MatDenseCUDAGetArray(A,&da);
421:   PetscMPIIntCast(A->cmap->n,&n);
422:   PetscMPIIntCast(A->rmap->n,&m);
423:   PetscMPIIntCast(a->lda,&lda);
424:   PetscInfo2(A,"LU factor %d x %d on backend\n",m,n);
425:   if (!dA->d_fact_ipiv) {
426:     ccer = cudaMalloc((void**)&dA->d_fact_ipiv,n*sizeof(*dA->d_fact_ipiv));CHKERRCUDA(ccer);
427:   }
428:   if (!dA->fact_lwork) {
429:     cerr = cusolverDnXgetrf_bufferSize(handle,m,n,da,lda,&dA->fact_lwork);CHKERRCUSOLVER(cerr);
430:     ccer = cudaMalloc((void**)&dA->d_fact_work,dA->fact_lwork*sizeof(*dA->d_fact_work));CHKERRCUDA(ccer);
431:   }
432:   if (!dA->d_fact_info) {
433:     ccer = cudaMalloc((void**)&dA->d_fact_info,sizeof(*dA->d_fact_info));CHKERRCUDA(ccer);
434:   }
435:   PetscLogGpuTimeBegin();
436:   cerr = cusolverDnXgetrf(handle,m,n,da,lda,dA->d_fact_work,dA->d_fact_ipiv,dA->d_fact_info);CHKERRCUSOLVER(cerr);
437:   WaitForGPU();CHKERRCUDA(ierr);
438:   PetscLogGpuTimeEnd();
439:   MatDenseCUDARestoreArray(A,&da);
440: #if defined(PETSC_USE_DEBUG)
441:   ccer = cudaMemcpy(&info, dA->d_fact_info, sizeof(int), cudaMemcpyDeviceToHost);CHKERRCUDA(ccer);
442:   if (info > 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_MAT_LU_ZRPVT,"Bad factorization: zero pivot in row %d",info-1);
443:   else if (info < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Wrong argument to cuSolver %d",-info);
444: #endif
445:   A->factortype = MAT_FACTOR_LU;
446:   PetscLogGpuFlops(2.0*n*n*m/3.0);

448:   A->ops->solve          = MatSolve_SeqDenseCUDA;
449:   A->ops->solvetranspose = MatSolveTranspose_SeqDenseCUDA;
450:   A->ops->matsolve       = MatMatSolve_SeqDenseCUDA;

452:   PetscFree(A->solvertype);
453:   PetscStrallocpy(MATSOLVERCUDA,&A->solvertype);
454:   return(0);
455: }

457: static PetscErrorCode MatCholeskyFactor_SeqDenseCUDA(Mat A,IS perm,const MatFactorInfo *factinfo)
458: {
459:   Mat_SeqDense       *a = (Mat_SeqDense*)A->data;
460:   Mat_SeqDenseCUDA   *dA = (Mat_SeqDenseCUDA*)A->spptr;
461:   PetscScalar        *da;
462:   int                n,lda;
463: #if defined(PETSC_USE_DEBUG)
464:   int                info;
465: #endif
466:   cusolverStatus_t   cerr;
467:   cusolverDnHandle_t handle;
468:   cudaError_t        ccer;
469:   PetscErrorCode     ierr;

472:   if (!A->rmap->n || !A->cmap->n) return(0);
473:   PetscCUSOLVERDnGetHandle(&handle);
474:   PetscMPIIntCast(A->rmap->n,&n);
475:   PetscInfo2(A,"Cholesky factor %d x %d on backend\n",n,n);
476:   if (A->spd) {
477:     MatDenseCUDAGetArray(A,&da);
478:     PetscMPIIntCast(a->lda,&lda);
479:     if (!dA->fact_lwork) {
480:       cerr = cusolverDnXpotrf_bufferSize(handle,CUBLAS_FILL_MODE_LOWER,n,da,lda,&dA->fact_lwork);CHKERRCUSOLVER(cerr);
481:       ccer = cudaMalloc((void**)&dA->d_fact_work,dA->fact_lwork*sizeof(*dA->d_fact_work));CHKERRCUDA(ccer);
482:     }
483:     if (!dA->d_fact_info) {
484:       ccer = cudaMalloc((void**)&dA->d_fact_info,sizeof(*dA->d_fact_info));CHKERRCUDA(ccer);
485:     }
486:     PetscLogGpuTimeBegin();
487:     cerr = cusolverDnXpotrf(handle,CUBLAS_FILL_MODE_LOWER,n,da,lda,dA->d_fact_work,dA->fact_lwork,dA->d_fact_info);CHKERRCUSOLVER(cerr);
488:     WaitForGPU();CHKERRCUDA(ierr);
489:     PetscLogGpuTimeEnd();

491:     MatDenseCUDARestoreArray(A,&da);
492: #if defined(PETSC_USE_DEBUG)
493:     ccer = cudaMemcpy(&info, dA->d_fact_info, sizeof(int), cudaMemcpyDeviceToHost);CHKERRCUDA(ccer);
494:     if (info > 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_MAT_CH_ZRPVT,"Bad factorization: zero pivot in row %d",info-1);
495:     else if (info < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Wrong argument to cuSolver %d",-info);
496: #endif
497:     A->factortype = MAT_FACTOR_CHOLESKY;
498:     PetscLogGpuFlops(1.0*n*n*n/3.0);
499:   } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"cusolverDnsytrs unavailable. Use MAT_FACTOR_LU");
500: #if 0
501:     /* at the time of writing this interface (cuda 10.0), cusolverDn does not implement *sytrs and *hetr* routines
502:        The code below should work, and it can be activated when *sytrs routines will be available */
503:     if (!dA->d_fact_ipiv) {
504:       ccer = cudaMalloc((void**)&dA->d_fact_ipiv,n*sizeof(*dA->d_fact_ipiv));CHKERRCUDA(ccer);
505:     }
506:     if (!dA->fact_lwork) {
507:       cerr = cusolverDnXsytrf_bufferSize(handle,n,da,lda,&dA->fact_lwork);CHKERRCUSOLVER(cerr);
508:       ccer = cudaMalloc((void**)&dA->d_fact_work,dA->fact_lwork*sizeof(*dA->d_fact_work));CHKERRCUDA(ccer);
509:     }
510:     if (!dA->d_fact_info) {
511:       ccer = cudaMalloc((void**)&dA->d_fact_info,sizeof(*dA->d_fact_info));CHKERRCUDA(ccer);
512:     }
513:     PetscLogGpuTimeBegin();
514:     cerr = cusolverDnXsytrf(handle,CUBLAS_FILL_MODE_LOWER,n,da,lda,dA->d_fact_ipiv,dA->d_fact_work,dA->fact_lwork,dA->d_fact_info);CHKERRCUSOLVER(cerr);
515:     PetscLogGpuTimeEnd();
516: #endif

518:   A->ops->solve          = MatSolve_SeqDenseCUDA;
519:   A->ops->solvetranspose = MatSolveTranspose_SeqDenseCUDA;
520:   A->ops->matsolve       = MatMatSolve_SeqDenseCUDA;

522:   PetscFree(A->solvertype);
523:   PetscStrallocpy(MATSOLVERCUDA,&A->solvertype);
524:   return(0);
525: }

527: /* GEMM kernel: C = op(A)*op(B), tA, tB flag transposition */
528: static PetscErrorCode MatMatMultNumeric_SeqDenseCUDA_SeqDenseCUDA_Private(Mat A,Mat B,Mat C,PetscBool tA, PetscBool tB)
529: {
530:   Mat_SeqDense      *a = (Mat_SeqDense*)A->data;
531:   Mat_SeqDense      *b = (Mat_SeqDense*)B->data;
532:   Mat_SeqDense      *c = (Mat_SeqDense*)C->data;
533:   const PetscScalar *da,*db;
534:   PetscScalar       *dc;
535:   PetscScalar       one=1.0,zero=0.0;
536:   int               m,n,k,alda,blda,clda;
537:   PetscErrorCode    ierr;
538:   cublasHandle_t    cublasv2handle;
539:   cublasStatus_t    berr;

542:   PetscMPIIntCast(C->rmap->n,&m);
543:   PetscMPIIntCast(C->cmap->n,&n);
544:   if (tA) {
545:     PetscMPIIntCast(A->rmap->n,&k);
546:   } else {
547:     PetscMPIIntCast(A->cmap->n,&k);
548:   }
549:   if (!m || !n || !k) return(0);
550:   PetscInfo3(C,"Matrix-Matrix product %d x %d x %d on backend\n",m,k,n);
551:   MatDenseCUDAGetArrayRead(A,&da);
552:   MatDenseCUDAGetArrayRead(B,&db);
553:   MatDenseCUDAGetArrayWrite(C,&dc);
554:   PetscMPIIntCast(a->lda,&alda);
555:   PetscMPIIntCast(b->lda,&blda);
556:   PetscMPIIntCast(c->lda,&clda);
557:   PetscCUBLASGetHandle(&cublasv2handle);
558:   PetscLogGpuTimeBegin();
559:   berr = cublasXgemm(cublasv2handle,tA ? CUBLAS_OP_T : CUBLAS_OP_N,tB ? CUBLAS_OP_T : CUBLAS_OP_N,
560:                      m,n,k,&one,da,alda,db,blda,&zero,dc,clda);CHKERRCUBLAS(berr);
561:   WaitForGPU();CHKERRCUDA(ierr);
562:   PetscLogGpuTimeEnd();
563:   PetscLogGpuFlops(1.0*m*n*k + 1.0*m*n*(k-1));
564:   MatDenseCUDARestoreArrayRead(A,&da);
565:   MatDenseCUDARestoreArrayRead(B,&db);
566:   MatDenseCUDARestoreArrayWrite(C,&dc);
567:   return(0);
568: }

570: PetscErrorCode MatTransposeMatMultNumeric_SeqDenseCUDA_SeqDenseCUDA(Mat A,Mat B,Mat C)
571: {

575:   MatMatMultNumeric_SeqDenseCUDA_SeqDenseCUDA_Private(A,B,C,PETSC_TRUE,PETSC_FALSE);
576:   return(0);
577: }

579: PetscErrorCode MatMatMultNumeric_SeqDenseCUDA_SeqDenseCUDA(Mat A,Mat B,Mat C)
580: {

584:   MatMatMultNumeric_SeqDenseCUDA_SeqDenseCUDA_Private(A,B,C,PETSC_FALSE,PETSC_FALSE);
585:   return(0);
586: }

588: PetscErrorCode MatMatTransposeMultNumeric_SeqDenseCUDA_SeqDenseCUDA(Mat A,Mat B,Mat C)
589: {

593:   MatMatMultNumeric_SeqDenseCUDA_SeqDenseCUDA_Private(A,B,C,PETSC_FALSE,PETSC_TRUE);
594:   return(0);
595: }

597: /* zz = op(A)*xx + yy
598:    if yy == NULL, only MatMult */
599: static PetscErrorCode MatMultAdd_SeqDenseCUDA_Private(Mat A,Vec xx,Vec yy,Vec zz,PetscBool trans)
600: {
601:   Mat_SeqDense      *mat = (Mat_SeqDense*)A->data;
602:   const PetscScalar *xarray,*da;
603:   PetscScalar       *zarray;
604:   PetscScalar       one=1.0,zero=0.0;
605:   int               m, n, lda; /* Use PetscMPIInt as it is typedef'ed to int */
606:   cublasHandle_t    cublasv2handle;
607:   cublasStatus_t    berr;
608:   PetscErrorCode    ierr;

611:   if (yy && yy != zz) { /* mult add */
612:     VecCopy_SeqCUDA(yy,zz);
613:   }
614:   if (!A->rmap->n || !A->cmap->n) {
615:     if (!yy) { /* mult only */
616:       VecSet_SeqCUDA(zz,0.0);
617:     }
618:     return(0);
619:   }
620:   PetscInfo2(A,"Matrix-vector product %d x %d on backend\n",A->rmap->n,A->cmap->n);
621:   PetscMPIIntCast(A->rmap->n,&m);
622:   PetscMPIIntCast(A->cmap->n,&n);
623:   PetscMPIIntCast(mat->lda,&lda);
624:   PetscCUBLASGetHandle(&cublasv2handle);
625:   MatDenseCUDAGetArrayRead(A,&da);
626:   VecCUDAGetArrayRead(xx,&xarray);
627:   VecCUDAGetArray(zz,&zarray);
628:   PetscLogGpuTimeBegin();
629:   berr = cublasXgemv(cublasv2handle,trans ? CUBLAS_OP_T : CUBLAS_OP_N,
630:                      m,n,&one,da,lda,xarray,1,(yy ? &one : &zero),zarray,1);CHKERRCUBLAS(berr);
631:   PetscLogGpuTimeEnd();
632:   PetscLogGpuFlops(2.0*A->rmap->n*A->cmap->n - (yy ? 0 : A->rmap->n));
633:   VecCUDARestoreArrayRead(xx,&xarray);
634:   VecCUDARestoreArray(zz,&zarray);
635:   MatDenseCUDARestoreArrayRead(A,&da);
636:   return(0);
637: }

639: PetscErrorCode MatMultAdd_SeqDenseCUDA(Mat A,Vec xx,Vec yy,Vec zz)
640: {

644:   MatMultAdd_SeqDenseCUDA_Private(A,xx,yy,zz,PETSC_FALSE);
645:   return(0);
646: }

648: PetscErrorCode MatMultTransposeAdd_SeqDenseCUDA(Mat A,Vec xx,Vec yy,Vec zz)
649: {

653:   MatMultAdd_SeqDenseCUDA_Private(A,xx,yy,zz,PETSC_TRUE);
654:   return(0);
655: }

657: PetscErrorCode MatMult_SeqDenseCUDA(Mat A,Vec xx,Vec yy)
658: {

662:   MatMultAdd_SeqDenseCUDA_Private(A,xx,NULL,yy,PETSC_FALSE);
663:   return(0);
664: }

666: PetscErrorCode MatMultTranspose_SeqDenseCUDA(Mat A,Vec xx,Vec yy)
667: {

671:   MatMultAdd_SeqDenseCUDA_Private(A,xx,NULL,yy,PETSC_TRUE);
672:   return(0);
673: }

675: PetscErrorCode MatDenseGetArrayRead_SeqDenseCUDA(Mat A,const PetscScalar *array[])
676: {
677:   Mat_SeqDense   *mat = (Mat_SeqDense*)A->data;

681:   MatSeqDenseCUDACopyFromGPU(A);
682:   *array = mat->v;
683:   return(0);
684: }

686: PetscErrorCode MatDenseGetArray_SeqDenseCUDA(Mat A,PetscScalar *array[])
687: {
688:   Mat_SeqDense   *mat = (Mat_SeqDense*)A->data;

692:   MatSeqDenseCUDACopyFromGPU(A);
693:   *array = mat->v;
694:   A->offloadmask = PETSC_OFFLOAD_CPU;
695:   return(0);
696: }

698: PetscErrorCode MatDenseRestoreArray_SeqDenseCUDA(Mat A,PetscScalar *array[])
699: {
701:   return(0);
702: }

704: PetscErrorCode MatAXPY_SeqDenseCUDA(Mat Y,PetscScalar alpha,Mat X,MatStructure str)
705: {
706:   Mat_SeqDense      *x = (Mat_SeqDense*)X->data;
707:   Mat_SeqDense      *y = (Mat_SeqDense*)Y->data;
708:   const PetscScalar *dx;
709:   PetscScalar       *dy;
710:   int               j,N,m,ldax,lday,one = 1;
711:   cublasHandle_t    cublasv2handle;
712:   cublasStatus_t    berr;
713:   PetscErrorCode    ierr;

716:   if (!X->rmap->n || !X->cmap->n) return(0);
717:   PetscCUBLASGetHandle(&cublasv2handle);
718:   MatDenseCUDAGetArrayRead(X,&dx);
719:   if (alpha != 0.0) {
720:     MatDenseCUDAGetArray(Y,&dy);
721:   } else {
722:     MatDenseCUDAGetArrayWrite(Y,&dy);
723:   }
724:   PetscMPIIntCast(X->rmap->n*X->cmap->n,&N);
725:   PetscMPIIntCast(X->rmap->n,&m);
726:   PetscMPIIntCast(x->lda,&ldax);
727:   PetscMPIIntCast(y->lda,&lday);
728:   PetscInfo2(Y,"Performing AXPY %d x %d on backend\n",Y->rmap->n,Y->cmap->n);
729:   PetscLogGpuTimeBegin();
730:   if (ldax>m || lday>m) {
731:     for (j=0; j<X->cmap->n; j++) {
732:       berr = cublasXaxpy(cublasv2handle,m,&alpha,dx+j*ldax,one,dy+j*lday,one);CHKERRCUBLAS(berr);
733:     }
734:   } else {
735:     berr = cublasXaxpy(cublasv2handle,N,&alpha,dx,one,dy,one);CHKERRCUBLAS(berr);
736:   }
737:   WaitForGPU();CHKERRCUDA(ierr);
738:   PetscLogGpuTimeEnd();
739:   PetscLogGpuFlops(PetscMax(2.*N-1,0));
740:   MatDenseCUDARestoreArrayRead(X,&dx);
741:   if (alpha != 0.0) {
742:     MatDenseCUDARestoreArray(Y,&dy);
743:   } else {
744:     MatDenseCUDARestoreArrayWrite(Y,&dy);
745:   }
746:   return(0);
747: }

749: static PetscErrorCode MatReset_SeqDenseCUDA(Mat A)
750: {
751:   Mat_SeqDenseCUDA *dA = (Mat_SeqDenseCUDA*)A->spptr;
752:   cudaError_t      cerr;
753:   PetscErrorCode   ierr;

756:   if (dA) {
757:     cerr = cudaFree(dA->d_v);CHKERRCUDA(cerr);
758:     cerr = cudaFree(dA->d_fact_ipiv);CHKERRCUDA(cerr);
759:     cerr = cudaFree(dA->d_fact_info);CHKERRCUDA(cerr);
760:     cerr = cudaFree(dA->d_fact_work);CHKERRCUDA(cerr);
761:     VecDestroy(&dA->workvec);
762:   }
763:   PetscFree(A->spptr);
764:   return(0);
765: }

767: PetscErrorCode MatDestroy_SeqDenseCUDA(Mat A)
768: {
769:   Mat_SeqDense   *a = (Mat_SeqDense*)A->data;

773:   /* prevent to copy back data if we own the data pointer */
774:   if (!a->user_alloc) { A->offloadmask = PETSC_OFFLOAD_CPU; }
775:   MatConvert_SeqDenseCUDA_SeqDense(A,MATSEQDENSE,MAT_INPLACE_MATRIX,&A);
776:   MatDestroy_SeqDense(A);
777:   return(0);
778: }

780: PetscErrorCode MatSeqDenseSetPreallocation_SeqDenseCUDA(Mat B,PetscScalar *data)
781: {
782:   Mat_SeqDense     *b;
783:   Mat_SeqDenseCUDA *dB;
784:   cudaError_t      cerr;
785:   PetscErrorCode   ierr;

788:   PetscLayoutSetUp(B->rmap);
789:   PetscLayoutSetUp(B->cmap);
790:   b       = (Mat_SeqDense*)B->data;
791:   b->Mmax = B->rmap->n;
792:   b->Nmax = B->cmap->n;
793:   if (b->lda <= 0 || b->changelda) b->lda = B->rmap->n;
794:   if (b->lda < B->rmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Invalid lda %D < %D",b->lda,B->rmap->n);

796:   PetscIntMultError(b->lda,b->Nmax,NULL);

798:   MatReset_SeqDenseCUDA(B);
799:   PetscNewLog(B,&dB);
800:   B->spptr = dB;
801:   cerr     = cudaMalloc((void**)&dB->d_v,b->lda*b->Nmax*sizeof(PetscScalar));CHKERRCUDA(cerr);

803:   if (!data) { /* petsc-allocated storage */
804:     if (!b->user_alloc) { PetscFree(b->v); }
805:     PetscCalloc1((size_t)b->lda*b->Nmax,&b->v);
806:     PetscLogObjectMemory((PetscObject)B,b->lda*b->Nmax*sizeof(PetscScalar));
807:     b->user_alloc       = PETSC_FALSE;
808:   } else { /* user-allocated storage */
809:     if (!b->user_alloc) { PetscFree(b->v); }
810:     b->v                = data;
811:     b->user_alloc       = PETSC_TRUE;
812:   }
813:   B->offloadmask = PETSC_OFFLOAD_CPU;
814:   B->preallocated     = PETSC_TRUE;
815:   B->assembled        = PETSC_TRUE;
816:   return(0);
817: }

819: PetscErrorCode MatDuplicate_SeqDenseCUDA(Mat A,MatDuplicateOption cpvalues,Mat *B)
820: {

824:   MatCreate(PetscObjectComm((PetscObject)A),B);
825:   MatSetSizes(*B,A->rmap->n,A->cmap->n,A->rmap->n,A->cmap->n);
826:   MatSetType(*B,((PetscObject)A)->type_name);
827:   MatDuplicateNoCreate_SeqDense(*B,A,cpvalues);
828:   if (cpvalues == MAT_COPY_VALUES && A->offloadmask != PETSC_OFFLOAD_CPU) {
829:     Mat_SeqDense      *a = (Mat_SeqDense*)A->data;
830:     const PetscScalar *da;
831:     PetscScalar       *db;
832:     cudaError_t       cerr;

834:     MatDenseCUDAGetArrayRead(A,&da);
835:     MatDenseCUDAGetArrayWrite(*B,&db);
836:     if (a->lda > A->rmap->n) {
837:       PetscInt j,m = A->rmap->n;

839:       for (j=0; j<A->cmap->n; j++) { /* it can be done better */
840:         cerr = cudaMemcpy(db+j*m,da+j*a->lda,m*sizeof(PetscScalar),cudaMemcpyDeviceToDevice);CHKERRCUDA(cerr);
841:       }
842:     } else {
843:       cerr = cudaMemcpy(db,da,a->lda*sizeof(PetscScalar)*A->cmap->n,cudaMemcpyDeviceToDevice);CHKERRCUDA(cerr);
844:     }
845:     MatDenseCUDARestoreArrayRead(A,&da);
846:     MatDenseCUDARestoreArrayWrite(*B,&db);
847:     (*B)->offloadmask = PETSC_OFFLOAD_BOTH;
848:   }
849:   return(0);
850: }

852: PETSC_INTERN PetscErrorCode MatGetFactor_seqdense_cuda(Mat A,MatFactorType ftype,Mat *fact)
853: {

857:   MatCreate(PetscObjectComm((PetscObject)A),fact);
858:   MatSetSizes(*fact,A->rmap->n,A->cmap->n,A->rmap->n,A->cmap->n);
859:   MatSetType(*fact,MATSEQDENSECUDA);
860:   if (ftype == MAT_FACTOR_LU) {
861:     (*fact)->ops->lufactorsymbolic = MatLUFactorSymbolic_SeqDense;
862:   } else {
863:     (*fact)->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_SeqDense;
864:   }
865:   (*fact)->factortype = ftype;

867:   PetscFree((*fact)->solvertype);
868:   PetscStrallocpy(MATSOLVERCUDA,&(*fact)->solvertype);
869:   return(0);
870: }

872: static PetscErrorCode MatPinToCPU_SeqDenseCUDA(Mat A,PetscBool flg)
873: {

877:   A->pinnedtocpu = flg;
878:   if (!flg) {
879:     PetscObjectComposeFunction((PetscObject)A,"MatSeqDenseSetPreallocation_C",MatSeqDenseSetPreallocation_SeqDenseCUDA);
880:     PetscObjectComposeFunction((PetscObject)A,"MatDenseGetArray_C",           MatDenseGetArray_SeqDenseCUDA);
881:     PetscObjectComposeFunction((PetscObject)A,"MatDenseGetArrayRead_C",       MatDenseGetArrayRead_SeqDenseCUDA);
882:     PetscObjectComposeFunction((PetscObject)A,"MatDenseRestoreArray_C",       MatDenseRestoreArray_SeqDenseCUDA);

884:     A->ops->duplicate               = MatDuplicate_SeqDenseCUDA;
885:     A->ops->mult                    = MatMult_SeqDenseCUDA;
886:     A->ops->multadd                 = MatMultAdd_SeqDenseCUDA;
887:     A->ops->multtranspose           = MatMultTranspose_SeqDenseCUDA;
888:     A->ops->multtransposeadd        = MatMultTransposeAdd_SeqDenseCUDA;
889:     A->ops->matmultnumeric          = MatMatMultNumeric_SeqDenseCUDA_SeqDenseCUDA;
890:     A->ops->mattransposemultnumeric = MatMatTransposeMultNumeric_SeqDenseCUDA_SeqDenseCUDA;
891:     A->ops->transposematmultnumeric = MatTransposeMatMultNumeric_SeqDenseCUDA_SeqDenseCUDA;
892:     A->ops->axpy                    = MatAXPY_SeqDenseCUDA;
893:     A->ops->choleskyfactor          = MatCholeskyFactor_SeqDenseCUDA;
894:     A->ops->lufactor                = MatLUFactor_SeqDenseCUDA;
895:   } else {
896:     /* make sure we have an up-to-date copy on the CPU */
897:     MatSeqDenseCUDACopyFromGPU(A);
898:     PetscObjectComposeFunction((PetscObject)A,"MatSeqDenseSetPreallocation_C",MatSeqDenseSetPreallocation_SeqDense);
899:     PetscObjectComposeFunction((PetscObject)A,"MatDenseGetArray_C",           MatDenseGetArray_SeqDense);
900:     PetscObjectComposeFunction((PetscObject)A,"MatDenseGetArrayRead_C",       MatDenseGetArray_SeqDense);
901:     PetscObjectComposeFunction((PetscObject)A,"MatDenseRestoreArray_C",       MatDenseRestoreArray_SeqDense);

903:     A->ops->duplicate               = MatDuplicate_SeqDense;
904:     A->ops->mult                    = MatMult_SeqDense;
905:     A->ops->multadd                 = MatMultAdd_SeqDense;
906:     A->ops->multtranspose           = MatMultTranspose_SeqDense;
907:     A->ops->multtransposeadd        = MatMultTransposeAdd_SeqDense;
908:     A->ops->matmultnumeric          = MatMatMultNumeric_SeqDense_SeqDense;
909:     A->ops->mattransposemultnumeric = MatMatTransposeMultNumeric_SeqDense_SeqDense;
910:     A->ops->transposematmultnumeric = MatTransposeMatMultNumeric_SeqDense_SeqDense;
911:     A->ops->axpy                    = MatAXPY_SeqDense;
912:     A->ops->choleskyfactor          = MatCholeskyFactor_SeqDense;
913:     A->ops->lufactor                = MatLUFactor_SeqDense;
914:  }
915:   return(0);
916: }

918: PetscErrorCode MatConvert_SeqDenseCUDA_SeqDense(Mat M,MatType type,MatReuse reuse,Mat *newmat)
919: {
920:   Mat              B;
921:   PetscErrorCode   ierr;

924:   if (reuse == MAT_REUSE_MATRIX || reuse == MAT_INITIAL_MATRIX) {
925:     /* TODO these cases should be optimized */
926:     MatConvert_Basic(M,type,reuse,newmat);
927:     return(0);
928:   }

930:   B    = *newmat;
931:   MatPinToCPU_SeqDenseCUDA(B,PETSC_TRUE);
932:   MatReset_SeqDenseCUDA(B);
933:   PetscFree(B->defaultvectype);
934:   PetscStrallocpy(VECSTANDARD,&B->defaultvectype);
935:   PetscObjectChangeTypeName((PetscObject)B,MATSEQDENSE);
936:   PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqdensecuda_seqdense_C",NULL);

938:   B->ops->pintocpu    = NULL;
939:   B->ops->destroy     = MatDestroy_SeqDense;
940:   B->offloadmask = PETSC_OFFLOAD_CPU;
941:   return(0);
942: }

944: PetscErrorCode MatConvert_SeqDense_SeqDenseCUDA(Mat M,MatType type,MatReuse reuse,Mat *newmat)
945: {
946:   Mat_SeqDenseCUDA *dB;
947:   Mat              B;
948:   PetscErrorCode   ierr;

951:   if (reuse == MAT_REUSE_MATRIX || reuse == MAT_INITIAL_MATRIX) {
952:     /* TODO these cases should be optimized */
953:     MatConvert_Basic(M,type,reuse,newmat);
954:     return(0);
955:   }

957:   B    = *newmat;
958:   PetscFree(B->defaultvectype);
959:   PetscStrallocpy(VECCUDA,&B->defaultvectype);
960:   PetscObjectChangeTypeName((PetscObject)B,MATSEQDENSECUDA);
961:   PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqdensecuda_seqdense_C",MatConvert_SeqDenseCUDA_SeqDense);

963:   MatReset_SeqDenseCUDA(B);
964:   PetscNewLog(B,&dB);
965:   B->spptr = dB;

967:   B->offloadmask = PETSC_OFFLOAD_UNALLOCATED;

969:   MatPinToCPU_SeqDenseCUDA(B,PETSC_FALSE);
970:   B->ops->pintocpu = MatPinToCPU_SeqDenseCUDA;
971:   B->ops->destroy  = MatDestroy_SeqDenseCUDA;
972:   return(0);
973: }

975: /*MC
976:    MATSEQDENSECUDA - MATSEQDENSECUDA = "seqdensecuda" - A matrix type to be used for sequential dense matrices on GPUs.

978:    Options Database Keys:
979: . -mat_type seqdensecuda - sets the matrix type to "seqdensecuda" during a call to MatSetFromOptions()

981:   Level: beginner

983: .seealso: MatCreateSeqDenseCuda()

985: M*/
986: PETSC_EXTERN PetscErrorCode MatCreate_SeqDenseCUDA(Mat B)
987: {

991:   MatCreate_SeqDense(B);
992:   MatConvert_SeqDense_SeqDenseCUDA(B,MATSEQDENSECUDA,MAT_INPLACE_MATRIX,&B);
993:   return(0);
994: }