Actual source code: densecuda.cu

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

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

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

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

392:   MatSolve_SeqDenseCUDA_Private(A,xx,yy,PETSC_FALSE);
393:   return(0);
394: }

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

401:   MatSolve_SeqDenseCUDA_Private(A,xx,yy,PETSC_TRUE);
402:   return(0);
403: }

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

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

450:   A->ops->solve          = MatSolve_SeqDenseCUDA;
451:   A->ops->solvetranspose = MatSolveTranspose_SeqDenseCUDA;
452:   A->ops->matsolve       = MatMatSolve_SeqDenseCUDA;

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

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

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

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

520:   A->ops->solve          = MatSolve_SeqDenseCUDA;
521:   A->ops->solvetranspose = MatSolveTranspose_SeqDenseCUDA;
522:   A->ops->matsolve       = MatMatSolve_SeqDenseCUDA;

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

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

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

573: PetscErrorCode MatTransposeMatMultNumeric_SeqDenseCUDA_SeqDenseCUDA(Mat A,Mat B,Mat C)
574: {

578:   MatMatMultNumeric_SeqDenseCUDA_SeqDenseCUDA_Private(A,B,C,PETSC_TRUE,PETSC_FALSE);
579:   return(0);
580: }

582: PetscErrorCode MatMatMultNumeric_SeqDenseCUDA_SeqDenseCUDA(Mat A,Mat B,Mat C)
583: {

587:   MatMatMultNumeric_SeqDenseCUDA_SeqDenseCUDA_Private(A,B,C,PETSC_FALSE,PETSC_FALSE);
588:   return(0);
589: }

591: PetscErrorCode MatMatTransposeMultNumeric_SeqDenseCUDA_SeqDenseCUDA(Mat A,Mat B,Mat C)
592: {

596:   MatMatMultNumeric_SeqDenseCUDA_SeqDenseCUDA_Private(A,B,C,PETSC_FALSE,PETSC_TRUE);
597:   return(0);
598: }

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

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

642: PetscErrorCode MatMultAdd_SeqDenseCUDA(Mat A,Vec xx,Vec yy,Vec zz)
643: {

647:   MatMultAdd_SeqDenseCUDA_Private(A,xx,yy,zz,PETSC_FALSE);
648:   return(0);
649: }

651: PetscErrorCode MatMultTransposeAdd_SeqDenseCUDA(Mat A,Vec xx,Vec yy,Vec zz)
652: {

656:   MatMultAdd_SeqDenseCUDA_Private(A,xx,yy,zz,PETSC_TRUE);
657:   return(0);
658: }

660: PetscErrorCode MatMult_SeqDenseCUDA(Mat A,Vec xx,Vec yy)
661: {

665:   MatMultAdd_SeqDenseCUDA_Private(A,xx,NULL,yy,PETSC_FALSE);
666:   return(0);
667: }

669: PetscErrorCode MatMultTranspose_SeqDenseCUDA(Mat A,Vec xx,Vec yy)
670: {

674:   MatMultAdd_SeqDenseCUDA_Private(A,xx,NULL,yy,PETSC_TRUE);
675:   return(0);
676: }

678: PetscErrorCode MatDenseGetArrayRead_SeqDenseCUDA(Mat A,const PetscScalar *array[])
679: {
680:   Mat_SeqDense   *mat = (Mat_SeqDense*)A->data;

684:   MatSeqDenseCUDACopyFromGPU(A);
685:   *array = mat->v;
686:   return(0);
687: }

689: PetscErrorCode MatDenseGetArray_SeqDenseCUDA(Mat A,PetscScalar *array[])
690: {
691:   Mat_SeqDense   *mat = (Mat_SeqDense*)A->data;

695:   MatSeqDenseCUDACopyFromGPU(A);
696:   *array = mat->v;
697:   A->offloadmask = PETSC_OFFLOAD_CPU;
698:   return(0);
699: }

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

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

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

753: static PetscErrorCode MatReset_SeqDenseCUDA(Mat A)
754: {
755:   Mat_SeqDenseCUDA *dA = (Mat_SeqDenseCUDA*)A->spptr;
756:   cudaError_t      cerr;
757:   PetscErrorCode   ierr;

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

771: PetscErrorCode MatDestroy_SeqDenseCUDA(Mat A)
772: {
773:   Mat_SeqDense   *a = (Mat_SeqDense*)A->data;

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

784: PetscErrorCode MatSeqDenseSetPreallocation_SeqDenseCUDA(Mat B,PetscScalar *data)
785: {
786:   Mat_SeqDense     *b;
787:   Mat_SeqDenseCUDA *dB;
788:   cudaError_t      cerr;
789:   PetscErrorCode   ierr;

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

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

802:   MatReset_SeqDenseCUDA(B);
803:   PetscNewLog(B,&dB);
804:   B->spptr = dB;
805:   cerr     = cudaMalloc((void**)&dB->d_v,b->lda*b->Nmax*sizeof(PetscScalar));CHKERRCUDA(cerr);

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

823: PetscErrorCode MatDuplicate_SeqDenseCUDA(Mat A,MatDuplicateOption cpvalues,Mat *B)
824: {

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

838:     MatDenseCUDAGetArrayRead(A,&da);
839:     MatDenseCUDAGetArrayWrite(*B,&db);
840:     if (a->lda > A->rmap->n) {
841:       PetscInt j,m = A->rmap->n;

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

856: PETSC_INTERN PetscErrorCode MatGetFactor_seqdense_cuda(Mat A,MatFactorType ftype,Mat *fact)
857: {

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

871:   PetscFree((*fact)->solvertype);
872:   PetscStrallocpy(MATSOLVERCUDA,&(*fact)->solvertype);
873:   return(0);
874: }

876: static PetscErrorCode MatPinToCPU_SeqDenseCUDA(Mat A,PetscBool flg)
877: {

881:   A->pinnedtocpu = flg;
882:   if (!flg) {
883:     PetscObjectComposeFunction((PetscObject)A,"MatSeqDenseSetPreallocation_C",MatSeqDenseSetPreallocation_SeqDenseCUDA);
884:     PetscObjectComposeFunction((PetscObject)A,"MatDenseGetArray_C",           MatDenseGetArray_SeqDenseCUDA);
885:     PetscObjectComposeFunction((PetscObject)A,"MatDenseGetArrayRead_C",       MatDenseGetArrayRead_SeqDenseCUDA);
886:     PetscObjectComposeFunction((PetscObject)A,"MatDenseRestoreArray_C",       MatDenseRestoreArray_SeqDenseCUDA);

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

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

922: PetscErrorCode MatConvert_SeqDenseCUDA_SeqDense(Mat M,MatType type,MatReuse reuse,Mat *newmat)
923: {
924:   Mat              B;
925:   PetscErrorCode   ierr;

928:   if (reuse == MAT_REUSE_MATRIX || reuse == MAT_INITIAL_MATRIX) {
929:     /* TODO these cases should be optimized */
930:     MatConvert_Basic(M,type,reuse,newmat);
931:     return(0);
932:   }

934:   B    = *newmat;
935:   MatPinToCPU_SeqDenseCUDA(B,PETSC_TRUE);
936:   MatReset_SeqDenseCUDA(B);
937:   PetscFree(B->defaultvectype);
938:   PetscStrallocpy(VECSTANDARD,&B->defaultvectype);
939:   PetscObjectChangeTypeName((PetscObject)B,MATSEQDENSE);
940:   PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqdensecuda_seqdense_C",NULL);

942:   B->ops->pintocpu    = NULL;
943:   B->ops->destroy     = MatDestroy_SeqDense;
944:   B->offloadmask = PETSC_OFFLOAD_CPU;
945:   return(0);
946: }

948: PetscErrorCode MatConvert_SeqDense_SeqDenseCUDA(Mat M,MatType type,MatReuse reuse,Mat *newmat)
949: {
950:   Mat_SeqDenseCUDA *dB;
951:   Mat              B;
952:   PetscErrorCode   ierr;

955:   if (reuse == MAT_REUSE_MATRIX || reuse == MAT_INITIAL_MATRIX) {
956:     /* TODO these cases should be optimized */
957:     MatConvert_Basic(M,type,reuse,newmat);
958:     return(0);
959:   }

961:   B    = *newmat;
962:   PetscFree(B->defaultvectype);
963:   PetscStrallocpy(VECCUDA,&B->defaultvectype);
964:   PetscObjectChangeTypeName((PetscObject)B,MATSEQDENSECUDA);
965:   PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqdensecuda_seqdense_C",MatConvert_SeqDenseCUDA_SeqDense);

967:   MatReset_SeqDenseCUDA(B);
968:   PetscNewLog(B,&dB);
969:   B->spptr = dB;

971:   B->offloadmask = PETSC_OFFLOAD_UNALLOCATED;

973:   MatPinToCPU_SeqDenseCUDA(B,PETSC_FALSE);
974:   B->ops->pintocpu = MatPinToCPU_SeqDenseCUDA;
975:   B->ops->destroy  = MatDestroy_SeqDenseCUDA;
976:   return(0);
977: }

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

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

985:   Level: beginner

987: .seealso: MatCreateSeqDenseCuda()

989: M*/
990: PETSC_EXTERN PetscErrorCode MatCreate_SeqDenseCUDA(Mat B)
991: {

995:   MatCreate_SeqDense(B);
996:   MatConvert_SeqDense_SeqDenseCUDA(B,MATSEQDENSECUDA,MAT_INPLACE_MATRIX,&B);
997:   return(0);
998: }