1 2 /* 3 Defines the basic matrix operations for sequential dense. 4 */ 5 6 #include <../src/mat/impls/dense/seq/dense.h> /*I "petscmat.h" I*/ 7 #include <petscblaslapack.h> 8 9 #include <../src/mat/impls/aij/seq/aij.h> 10 11 PetscErrorCode MatSeqDenseSymmetrize_Private(Mat A, PetscBool hermitian) 12 { 13 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 14 PetscInt j, k, n = A->rmap->n; 15 PetscScalar *v; 16 PetscErrorCode ierr; 17 18 PetscFunctionBegin; 19 if (A->rmap->n != A->cmap->n) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Cannot symmetrize a rectangular matrix"); 20 ierr = MatDenseGetArray(A,&v);CHKERRQ(ierr); 21 if (!hermitian) { 22 for (k=0;k<n;k++) { 23 for (j=k;j<n;j++) { 24 v[j*mat->lda + k] = v[k*mat->lda + j]; 25 } 26 } 27 } else { 28 for (k=0;k<n;k++) { 29 for (j=k;j<n;j++) { 30 v[j*mat->lda + k] = PetscConj(v[k*mat->lda + j]); 31 } 32 } 33 } 34 ierr = MatDenseRestoreArray(A,&v);CHKERRQ(ierr); 35 PetscFunctionReturn(0); 36 } 37 38 PETSC_EXTERN PetscErrorCode MatSeqDenseInvertFactors_Private(Mat A) 39 { 40 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 41 PetscErrorCode ierr; 42 PetscBLASInt info,n; 43 44 PetscFunctionBegin; 45 if (!A->rmap->n || !A->cmap->n) PetscFunctionReturn(0); 46 ierr = PetscBLASIntCast(A->cmap->n,&n);CHKERRQ(ierr); 47 if (A->factortype == MAT_FACTOR_LU) { 48 if (!mat->pivots) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Pivots not present"); 49 if (!mat->fwork) { 50 mat->lfwork = n; 51 ierr = PetscMalloc1(mat->lfwork,&mat->fwork);CHKERRQ(ierr); 52 ierr = PetscLogObjectMemory((PetscObject)A,mat->lfwork*sizeof(PetscBLASInt));CHKERRQ(ierr); 53 } 54 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr); 55 PetscStackCallBLAS("LAPACKgetri",LAPACKgetri_(&n,mat->v,&mat->lda,mat->pivots,mat->fwork,&mat->lfwork,&info)); 56 ierr = PetscFPTrapPop();CHKERRQ(ierr); 57 ierr = PetscLogFlops((1.0*A->cmap->n*A->cmap->n*A->cmap->n)/3.0);CHKERRQ(ierr); 58 } else if (A->factortype == MAT_FACTOR_CHOLESKY) { 59 if (A->spd) { 60 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr); 61 PetscStackCallBLAS("LAPACKpotri",LAPACKpotri_("L",&n,mat->v,&mat->lda,&info)); 62 ierr = PetscFPTrapPop();CHKERRQ(ierr); 63 ierr = MatSeqDenseSymmetrize_Private(A,PETSC_TRUE);CHKERRQ(ierr); 64 #if defined(PETSC_USE_COMPLEX) 65 } else if (A->hermitian) { 66 if (!mat->pivots) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Pivots not present"); 67 if (!mat->fwork) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Fwork not present"); 68 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr); 69 PetscStackCallBLAS("LAPACKhetri",LAPACKhetri_("L",&n,mat->v,&mat->lda,mat->pivots,mat->fwork,&info)); 70 ierr = PetscFPTrapPop();CHKERRQ(ierr); 71 ierr = MatSeqDenseSymmetrize_Private(A,PETSC_TRUE);CHKERRQ(ierr); 72 #endif 73 } else { /* symmetric case */ 74 if (!mat->pivots) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Pivots not present"); 75 if (!mat->fwork) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Fwork not present"); 76 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr); 77 PetscStackCallBLAS("LAPACKsytri",LAPACKsytri_("L",&n,mat->v,&mat->lda,mat->pivots,mat->fwork,&info)); 78 ierr = PetscFPTrapPop();CHKERRQ(ierr); 79 ierr = MatSeqDenseSymmetrize_Private(A,PETSC_FALSE);CHKERRQ(ierr); 80 } 81 if (info) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_MAT_CH_ZRPVT,"Bad Inversion: zero pivot in row %D",(PetscInt)info-1); 82 ierr = PetscLogFlops((1.0*A->cmap->n*A->cmap->n*A->cmap->n)/3.0);CHKERRQ(ierr); 83 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix must be factored to solve"); 84 85 A->ops->solve = NULL; 86 A->ops->matsolve = NULL; 87 A->ops->solvetranspose = NULL; 88 A->ops->matsolvetranspose = NULL; 89 A->ops->solveadd = NULL; 90 A->ops->solvetransposeadd = NULL; 91 A->factortype = MAT_FACTOR_NONE; 92 ierr = PetscFree(A->solvertype);CHKERRQ(ierr); 93 PetscFunctionReturn(0); 94 } 95 96 PetscErrorCode MatZeroRowsColumns_SeqDense(Mat A,PetscInt N,const PetscInt rows[],PetscScalar diag,Vec x,Vec b) 97 { 98 PetscErrorCode ierr; 99 Mat_SeqDense *l = (Mat_SeqDense*)A->data; 100 PetscInt m = l->lda, n = A->cmap->n,r = A->rmap->n, i,j; 101 PetscScalar *slot,*bb,*v; 102 const PetscScalar *xx; 103 104 PetscFunctionBegin; 105 if (PetscDefined(USE_DEBUG)) { 106 for (i=0; i<N; i++) { 107 if (rows[i] < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative row requested to be zeroed"); 108 if (rows[i] >= A->rmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row %D requested to be zeroed greater than or equal number of rows %D",rows[i],A->rmap->n); 109 if (rows[i] >= A->cmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Col %D requested to be zeroed greater than or equal number of cols %D",rows[i],A->cmap->n); 110 } 111 } 112 if (!N) PetscFunctionReturn(0); 113 114 /* fix right hand side if needed */ 115 if (x && b) { 116 Vec xt; 117 118 if (A->rmap->n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only coded for square matrices"); 119 ierr = VecDuplicate(x,&xt);CHKERRQ(ierr); 120 ierr = VecCopy(x,xt);CHKERRQ(ierr); 121 ierr = VecScale(xt,-1.0);CHKERRQ(ierr); 122 ierr = MatMultAdd(A,xt,b,b);CHKERRQ(ierr); 123 ierr = VecDestroy(&xt);CHKERRQ(ierr); 124 ierr = VecGetArrayRead(x,&xx);CHKERRQ(ierr); 125 ierr = VecGetArray(b,&bb);CHKERRQ(ierr); 126 for (i=0; i<N; i++) bb[rows[i]] = diag*xx[rows[i]]; 127 ierr = VecRestoreArrayRead(x,&xx);CHKERRQ(ierr); 128 ierr = VecRestoreArray(b,&bb);CHKERRQ(ierr); 129 } 130 131 ierr = MatDenseGetArray(A,&v);CHKERRQ(ierr); 132 for (i=0; i<N; i++) { 133 slot = v + rows[i]*m; 134 ierr = PetscArrayzero(slot,r);CHKERRQ(ierr); 135 } 136 for (i=0; i<N; i++) { 137 slot = v + rows[i]; 138 for (j=0; j<n; j++) { *slot = 0.0; slot += m;} 139 } 140 if (diag != 0.0) { 141 if (A->rmap->n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only coded for square matrices"); 142 for (i=0; i<N; i++) { 143 slot = v + (m+1)*rows[i]; 144 *slot = diag; 145 } 146 } 147 ierr = MatDenseRestoreArray(A,&v);CHKERRQ(ierr); 148 PetscFunctionReturn(0); 149 } 150 151 PetscErrorCode MatPtAPNumeric_SeqDense_SeqDense(Mat A,Mat P,Mat C) 152 { 153 Mat_SeqDense *c = (Mat_SeqDense*)(C->data); 154 PetscErrorCode ierr; 155 156 PetscFunctionBegin; 157 if (c->ptapwork) { 158 ierr = (*C->ops->matmultnumeric)(A,P,c->ptapwork);CHKERRQ(ierr); 159 ierr = (*C->ops->transposematmultnumeric)(P,c->ptapwork,C);CHKERRQ(ierr); 160 } else SETERRQ(PetscObjectComm((PetscObject)C),PETSC_ERR_SUP,"Must call MatPtAPSymbolic_SeqDense_SeqDense() first"); 161 PetscFunctionReturn(0); 162 } 163 164 PetscErrorCode MatPtAPSymbolic_SeqDense_SeqDense(Mat A,Mat P,PetscReal fill,Mat C) 165 { 166 Mat_SeqDense *c; 167 PetscBool flg; 168 PetscErrorCode ierr; 169 170 PetscFunctionBegin; 171 ierr = PetscObjectTypeCompare((PetscObject)P,((PetscObject)A)->type_name,&flg);CHKERRQ(ierr); 172 ierr = MatSetSizes(C,P->cmap->n,P->cmap->n,P->cmap->N,P->cmap->N);CHKERRQ(ierr); 173 ierr = MatSetType(C,flg ? ((PetscObject)A)->type_name : MATDENSE);CHKERRQ(ierr); 174 ierr = MatSeqDenseSetPreallocation(C,NULL);CHKERRQ(ierr); 175 c = (Mat_SeqDense*)C->data; 176 ierr = MatCreate(PetscObjectComm((PetscObject)A),&c->ptapwork);CHKERRQ(ierr); 177 ierr = MatSetSizes(c->ptapwork,A->rmap->n,P->cmap->n,A->rmap->N,P->cmap->N);CHKERRQ(ierr); 178 ierr = MatSetType(c->ptapwork,flg ? ((PetscObject)A)->type_name : MATDENSE);CHKERRQ(ierr); 179 ierr = MatSeqDenseSetPreallocation(c->ptapwork,NULL);CHKERRQ(ierr); 180 PetscFunctionReturn(0); 181 } 182 183 PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_SeqDense(Mat A,MatType newtype,MatReuse reuse,Mat *newmat) 184 { 185 Mat B = NULL; 186 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 187 Mat_SeqDense *b; 188 PetscErrorCode ierr; 189 PetscInt *ai=a->i,*aj=a->j,m=A->rmap->N,n=A->cmap->N,i; 190 MatScalar *av=a->a; 191 PetscBool isseqdense; 192 193 PetscFunctionBegin; 194 if (reuse == MAT_REUSE_MATRIX) { 195 ierr = PetscObjectTypeCompare((PetscObject)*newmat,MATSEQDENSE,&isseqdense);CHKERRQ(ierr); 196 if (!isseqdense) SETERRQ1(PetscObjectComm((PetscObject)*newmat),PETSC_ERR_USER,"Cannot reuse matrix of type %s",((PetscObject)(*newmat))->type_name); 197 } 198 if (reuse != MAT_REUSE_MATRIX) { 199 ierr = MatCreate(PetscObjectComm((PetscObject)A),&B);CHKERRQ(ierr); 200 ierr = MatSetSizes(B,m,n,m,n);CHKERRQ(ierr); 201 ierr = MatSetType(B,MATSEQDENSE);CHKERRQ(ierr); 202 ierr = MatSeqDenseSetPreallocation(B,NULL);CHKERRQ(ierr); 203 b = (Mat_SeqDense*)(B->data); 204 } else { 205 b = (Mat_SeqDense*)((*newmat)->data); 206 ierr = PetscArrayzero(b->v,m*n);CHKERRQ(ierr); 207 } 208 for (i=0; i<m; i++) { 209 PetscInt j; 210 for (j=0;j<ai[1]-ai[0];j++) { 211 b->v[*aj*m+i] = *av; 212 aj++; 213 av++; 214 } 215 ai++; 216 } 217 218 if (reuse == MAT_INPLACE_MATRIX) { 219 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 220 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 221 ierr = MatHeaderReplace(A,&B);CHKERRQ(ierr); 222 } else { 223 if (B) *newmat = B; 224 ierr = MatAssemblyBegin(*newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 225 ierr = MatAssemblyEnd(*newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 226 } 227 PetscFunctionReturn(0); 228 } 229 230 PETSC_INTERN PetscErrorCode MatConvert_SeqDense_SeqAIJ(Mat A, MatType newtype,MatReuse reuse,Mat *newmat) 231 { 232 Mat B; 233 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 234 PetscErrorCode ierr; 235 PetscInt i, j; 236 PetscInt *rows, *nnz; 237 MatScalar *aa = a->v, *vals; 238 239 PetscFunctionBegin; 240 ierr = MatCreate(PetscObjectComm((PetscObject)A),&B);CHKERRQ(ierr); 241 ierr = MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);CHKERRQ(ierr); 242 ierr = MatSetType(B,MATSEQAIJ);CHKERRQ(ierr); 243 ierr = PetscCalloc3(A->rmap->n,&rows,A->rmap->n,&nnz,A->rmap->n,&vals);CHKERRQ(ierr); 244 for (j=0; j<A->cmap->n; j++) { 245 for (i=0; i<A->rmap->n; i++) if (aa[i] != 0.0 || i == j) ++nnz[i]; 246 aa += a->lda; 247 } 248 ierr = MatSeqAIJSetPreallocation(B,PETSC_DETERMINE,nnz);CHKERRQ(ierr); 249 aa = a->v; 250 for (j=0; j<A->cmap->n; j++) { 251 PetscInt numRows = 0; 252 for (i=0; i<A->rmap->n; i++) if (aa[i] != 0.0 || i == j) {rows[numRows] = i; vals[numRows++] = aa[i];} 253 ierr = MatSetValues(B,numRows,rows,1,&j,vals,INSERT_VALUES);CHKERRQ(ierr); 254 aa += a->lda; 255 } 256 ierr = PetscFree3(rows,nnz,vals);CHKERRQ(ierr); 257 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 258 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 259 260 if (reuse == MAT_INPLACE_MATRIX) { 261 ierr = MatHeaderReplace(A,&B);CHKERRQ(ierr); 262 } else { 263 *newmat = B; 264 } 265 PetscFunctionReturn(0); 266 } 267 268 PetscErrorCode MatAXPY_SeqDense(Mat Y,PetscScalar alpha,Mat X,MatStructure str) 269 { 270 Mat_SeqDense *x = (Mat_SeqDense*)X->data,*y = (Mat_SeqDense*)Y->data; 271 const PetscScalar *xv; 272 PetscScalar *yv; 273 PetscBLASInt N,m,ldax,lday,one = 1; 274 PetscErrorCode ierr; 275 276 PetscFunctionBegin; 277 ierr = MatDenseGetArrayRead(X,&xv);CHKERRQ(ierr); 278 ierr = MatDenseGetArray(Y,&yv);CHKERRQ(ierr); 279 ierr = PetscBLASIntCast(X->rmap->n*X->cmap->n,&N);CHKERRQ(ierr); 280 ierr = PetscBLASIntCast(X->rmap->n,&m);CHKERRQ(ierr); 281 ierr = PetscBLASIntCast(x->lda,&ldax);CHKERRQ(ierr); 282 ierr = PetscBLASIntCast(y->lda,&lday);CHKERRQ(ierr); 283 if (ldax>m || lday>m) { 284 PetscInt j; 285 286 for (j=0; j<X->cmap->n; j++) { 287 PetscStackCallBLAS("BLASaxpy",BLASaxpy_(&m,&alpha,xv+j*ldax,&one,yv+j*lday,&one)); 288 } 289 } else { 290 PetscStackCallBLAS("BLASaxpy",BLASaxpy_(&N,&alpha,xv,&one,yv,&one)); 291 } 292 ierr = MatDenseRestoreArrayRead(X,&xv);CHKERRQ(ierr); 293 ierr = MatDenseRestoreArray(Y,&yv);CHKERRQ(ierr); 294 ierr = PetscLogFlops(PetscMax(2*N-1,0));CHKERRQ(ierr); 295 PetscFunctionReturn(0); 296 } 297 298 static PetscErrorCode MatGetInfo_SeqDense(Mat A,MatInfoType flag,MatInfo *info) 299 { 300 PetscLogDouble N = A->rmap->n*A->cmap->n; 301 302 PetscFunctionBegin; 303 info->block_size = 1.0; 304 info->nz_allocated = N; 305 info->nz_used = N; 306 info->nz_unneeded = 0; 307 info->assemblies = A->num_ass; 308 info->mallocs = 0; 309 info->memory = ((PetscObject)A)->mem; 310 info->fill_ratio_given = 0; 311 info->fill_ratio_needed = 0; 312 info->factor_mallocs = 0; 313 PetscFunctionReturn(0); 314 } 315 316 PetscErrorCode MatScale_SeqDense(Mat A,PetscScalar alpha) 317 { 318 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 319 PetscScalar *v; 320 PetscErrorCode ierr; 321 PetscBLASInt one = 1,j,nz,lda; 322 323 PetscFunctionBegin; 324 ierr = MatDenseGetArray(A,&v);CHKERRQ(ierr); 325 ierr = PetscBLASIntCast(a->lda,&lda);CHKERRQ(ierr); 326 if (lda>A->rmap->n) { 327 ierr = PetscBLASIntCast(A->rmap->n,&nz);CHKERRQ(ierr); 328 for (j=0; j<A->cmap->n; j++) { 329 PetscStackCallBLAS("BLASscal",BLASscal_(&nz,&alpha,v+j*lda,&one)); 330 } 331 } else { 332 ierr = PetscBLASIntCast(A->rmap->n*A->cmap->n,&nz);CHKERRQ(ierr); 333 PetscStackCallBLAS("BLASscal",BLASscal_(&nz,&alpha,v,&one)); 334 } 335 ierr = PetscLogFlops(nz);CHKERRQ(ierr); 336 ierr = MatDenseRestoreArray(A,&v);CHKERRQ(ierr); 337 PetscFunctionReturn(0); 338 } 339 340 static PetscErrorCode MatIsHermitian_SeqDense(Mat A,PetscReal rtol,PetscBool *fl) 341 { 342 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 343 PetscInt i,j,m = A->rmap->n,N = a->lda; 344 const PetscScalar *v; 345 PetscErrorCode ierr; 346 347 PetscFunctionBegin; 348 *fl = PETSC_FALSE; 349 if (A->rmap->n != A->cmap->n) PetscFunctionReturn(0); 350 ierr = MatDenseGetArrayRead(A,&v);CHKERRQ(ierr); 351 for (i=0; i<m; i++) { 352 for (j=i; j<m; j++) { 353 if (PetscAbsScalar(v[i+j*N] - PetscConj(v[j+i*N])) > rtol) { 354 goto restore; 355 } 356 } 357 } 358 *fl = PETSC_TRUE; 359 restore: 360 ierr = MatDenseRestoreArrayRead(A,&v);CHKERRQ(ierr); 361 PetscFunctionReturn(0); 362 } 363 364 static PetscErrorCode MatIsSymmetric_SeqDense(Mat A,PetscReal rtol,PetscBool *fl) 365 { 366 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 367 PetscInt i,j,m = A->rmap->n,N = a->lda; 368 const PetscScalar *v; 369 PetscErrorCode ierr; 370 371 PetscFunctionBegin; 372 *fl = PETSC_FALSE; 373 if (A->rmap->n != A->cmap->n) PetscFunctionReturn(0); 374 ierr = MatDenseGetArrayRead(A,&v);CHKERRQ(ierr); 375 for (i=0; i<m; i++) { 376 for (j=i; j<m; j++) { 377 if (PetscAbsScalar(v[i+j*N] - v[j+i*N]) > rtol) { 378 goto restore; 379 } 380 } 381 } 382 *fl = PETSC_TRUE; 383 restore: 384 ierr = MatDenseRestoreArrayRead(A,&v);CHKERRQ(ierr); 385 PetscFunctionReturn(0); 386 } 387 388 PetscErrorCode MatDuplicateNoCreate_SeqDense(Mat newi,Mat A,MatDuplicateOption cpvalues) 389 { 390 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 391 PetscErrorCode ierr; 392 PetscInt lda = (PetscInt)mat->lda,j,m; 393 394 PetscFunctionBegin; 395 ierr = PetscLayoutReference(A->rmap,&newi->rmap);CHKERRQ(ierr); 396 ierr = PetscLayoutReference(A->cmap,&newi->cmap);CHKERRQ(ierr); 397 ierr = MatSeqDenseSetPreallocation(newi,NULL);CHKERRQ(ierr); 398 if (cpvalues == MAT_COPY_VALUES) { 399 const PetscScalar *av; 400 PetscScalar *v; 401 402 ierr = MatDenseGetArrayRead(A,&av);CHKERRQ(ierr); 403 ierr = MatDenseGetArray(newi,&v);CHKERRQ(ierr); 404 if (lda>A->rmap->n) { 405 m = A->rmap->n; 406 for (j=0; j<A->cmap->n; j++) { 407 ierr = PetscArraycpy(v+j*m,av+j*lda,m);CHKERRQ(ierr); 408 } 409 } else { 410 ierr = PetscArraycpy(v,av,A->rmap->n*A->cmap->n);CHKERRQ(ierr); 411 } 412 ierr = MatDenseRestoreArray(newi,&v);CHKERRQ(ierr); 413 ierr = MatDenseRestoreArrayRead(A,&av);CHKERRQ(ierr); 414 } 415 PetscFunctionReturn(0); 416 } 417 418 PetscErrorCode MatDuplicate_SeqDense(Mat A,MatDuplicateOption cpvalues,Mat *newmat) 419 { 420 PetscErrorCode ierr; 421 422 PetscFunctionBegin; 423 ierr = MatCreate(PetscObjectComm((PetscObject)A),newmat);CHKERRQ(ierr); 424 ierr = MatSetSizes(*newmat,A->rmap->n,A->cmap->n,A->rmap->n,A->cmap->n);CHKERRQ(ierr); 425 ierr = MatSetType(*newmat,((PetscObject)A)->type_name);CHKERRQ(ierr); 426 ierr = MatDuplicateNoCreate_SeqDense(*newmat,A,cpvalues);CHKERRQ(ierr); 427 PetscFunctionReturn(0); 428 } 429 430 static PetscErrorCode MatLUFactorNumeric_SeqDense(Mat fact,Mat A,const MatFactorInfo *info_dummy) 431 { 432 MatFactorInfo info; 433 PetscErrorCode ierr; 434 435 PetscFunctionBegin; 436 ierr = MatDuplicateNoCreate_SeqDense(fact,A,MAT_COPY_VALUES);CHKERRQ(ierr); 437 ierr = (*fact->ops->lufactor)(fact,0,0,&info);CHKERRQ(ierr); 438 PetscFunctionReturn(0); 439 } 440 441 static PetscErrorCode MatSolve_SeqDense(Mat A,Vec xx,Vec yy) 442 { 443 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 444 PetscErrorCode ierr; 445 const PetscScalar *x; 446 PetscScalar *y; 447 PetscBLASInt one = 1,info,m; 448 449 PetscFunctionBegin; 450 ierr = PetscBLASIntCast(A->rmap->n,&m);CHKERRQ(ierr); 451 ierr = VecGetArrayRead(xx,&x);CHKERRQ(ierr); 452 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 453 ierr = PetscArraycpy(y,x,A->rmap->n);CHKERRQ(ierr); 454 if (A->factortype == MAT_FACTOR_LU) { 455 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr); 456 PetscStackCallBLAS("LAPACKgetrs",LAPACKgetrs_("N",&m,&one,mat->v,&mat->lda,mat->pivots,y,&m,&info)); 457 ierr = PetscFPTrapPop();CHKERRQ(ierr); 458 if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"GETRS - Bad solve"); 459 } else if (A->factortype == MAT_FACTOR_CHOLESKY) { 460 if (A->spd) { 461 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr); 462 PetscStackCallBLAS("LAPACKpotrs",LAPACKpotrs_("L",&m,&one,mat->v,&mat->lda,y,&m,&info)); 463 ierr = PetscFPTrapPop();CHKERRQ(ierr); 464 if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"POTRS Bad solve"); 465 #if defined(PETSC_USE_COMPLEX) 466 } else if (A->hermitian) { 467 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr); 468 PetscStackCallBLAS("LAPACKhetrs",LAPACKhetrs_("L",&m,&one,mat->v,&mat->lda,mat->pivots,y,&m,&info)); 469 ierr = PetscFPTrapPop();CHKERRQ(ierr); 470 if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"HETRS Bad solve"); 471 #endif 472 } else { /* symmetric case */ 473 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr); 474 PetscStackCallBLAS("LAPACKsytrs",LAPACKsytrs_("L",&m,&one,mat->v,&mat->lda,mat->pivots,y,&m,&info)); 475 ierr = PetscFPTrapPop();CHKERRQ(ierr); 476 if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"SYTRS Bad solve"); 477 } 478 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix must be factored to solve"); 479 ierr = VecRestoreArrayRead(xx,&x);CHKERRQ(ierr); 480 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 481 ierr = PetscLogFlops(2.0*A->cmap->n*A->cmap->n - A->cmap->n);CHKERRQ(ierr); 482 PetscFunctionReturn(0); 483 } 484 485 static PetscErrorCode MatMatSolve_SeqDense(Mat A,Mat B,Mat X) 486 { 487 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 488 PetscErrorCode ierr; 489 const PetscScalar *b; 490 PetscScalar *x; 491 PetscInt n; 492 PetscBLASInt nrhs,info,m; 493 494 PetscFunctionBegin; 495 ierr = PetscBLASIntCast(A->rmap->n,&m);CHKERRQ(ierr); 496 ierr = MatGetSize(B,NULL,&n);CHKERRQ(ierr); 497 ierr = PetscBLASIntCast(n,&nrhs);CHKERRQ(ierr); 498 ierr = MatDenseGetArrayRead(B,&b);CHKERRQ(ierr); 499 ierr = MatDenseGetArray(X,&x);CHKERRQ(ierr); 500 501 ierr = PetscArraycpy(x,b,m*nrhs);CHKERRQ(ierr); 502 503 if (A->factortype == MAT_FACTOR_LU) { 504 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr); 505 PetscStackCallBLAS("LAPACKgetrs",LAPACKgetrs_("N",&m,&nrhs,mat->v,&mat->lda,mat->pivots,x,&m,&info)); 506 ierr = PetscFPTrapPop();CHKERRQ(ierr); 507 if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"GETRS - Bad solve"); 508 } else if (A->factortype == MAT_FACTOR_CHOLESKY) { 509 if (A->spd) { 510 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr); 511 PetscStackCallBLAS("LAPACKpotrs",LAPACKpotrs_("L",&m,&nrhs,mat->v,&mat->lda,x,&m,&info)); 512 ierr = PetscFPTrapPop();CHKERRQ(ierr); 513 if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"POTRS Bad solve"); 514 #if defined(PETSC_USE_COMPLEX) 515 } else if (A->hermitian) { 516 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr); 517 PetscStackCallBLAS("LAPACKhetrs",LAPACKhetrs_("L",&m,&nrhs,mat->v,&mat->lda,mat->pivots,x,&m,&info)); 518 ierr = PetscFPTrapPop();CHKERRQ(ierr); 519 if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"HETRS Bad solve"); 520 #endif 521 } else { /* symmetric case */ 522 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr); 523 PetscStackCallBLAS("LAPACKsytrs",LAPACKsytrs_("L",&m,&nrhs,mat->v,&mat->lda,mat->pivots,x,&m,&info)); 524 ierr = PetscFPTrapPop();CHKERRQ(ierr); 525 if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"SYTRS Bad solve"); 526 } 527 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix must be factored to solve"); 528 529 ierr = MatDenseRestoreArrayRead(B,&b);CHKERRQ(ierr); 530 ierr = MatDenseRestoreArray(X,&x);CHKERRQ(ierr); 531 ierr = PetscLogFlops(nrhs*(2.0*m*m - m));CHKERRQ(ierr); 532 PetscFunctionReturn(0); 533 } 534 535 static PetscErrorCode MatConjugate_SeqDense(Mat); 536 537 static PetscErrorCode MatSolveTranspose_SeqDense(Mat A,Vec xx,Vec yy) 538 { 539 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 540 PetscErrorCode ierr; 541 const PetscScalar *x; 542 PetscScalar *y; 543 PetscBLASInt one = 1,info,m; 544 545 PetscFunctionBegin; 546 ierr = PetscBLASIntCast(A->rmap->n,&m);CHKERRQ(ierr); 547 ierr = VecGetArrayRead(xx,&x);CHKERRQ(ierr); 548 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 549 ierr = PetscArraycpy(y,x,A->rmap->n);CHKERRQ(ierr); 550 if (A->factortype == MAT_FACTOR_LU) { 551 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr); 552 PetscStackCallBLAS("LAPACKgetrs",LAPACKgetrs_("T",&m,&one,mat->v,&mat->lda,mat->pivots,y,&m,&info)); 553 ierr = PetscFPTrapPop();CHKERRQ(ierr); 554 if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"POTRS - Bad solve"); 555 } else if (A->factortype == MAT_FACTOR_CHOLESKY) { 556 if (A->spd) { 557 #if defined(PETSC_USE_COMPLEX) 558 ierr = MatConjugate_SeqDense(A);CHKERRQ(ierr); 559 #endif 560 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr); 561 PetscStackCallBLAS("LAPACKpotrs",LAPACKpotrs_("L",&m,&one,mat->v,&mat->lda,y,&m,&info)); 562 ierr = PetscFPTrapPop();CHKERRQ(ierr); 563 #if defined(PETSC_USE_COMPLEX) 564 ierr = MatConjugate_SeqDense(A);CHKERRQ(ierr); 565 #endif 566 if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"POTRS Bad solve"); 567 #if defined(PETSC_USE_COMPLEX) 568 } else if (A->hermitian) { 569 ierr = MatConjugate_SeqDense(A);CHKERRQ(ierr); 570 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr); 571 PetscStackCallBLAS("LAPACKhetrs",LAPACKhetrs_("L",&m,&one,mat->v,&mat->lda,mat->pivots,y,&m,&info)); 572 ierr = PetscFPTrapPop();CHKERRQ(ierr); 573 ierr = MatConjugate_SeqDense(A);CHKERRQ(ierr); 574 #endif 575 } else { /* symmetric case */ 576 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr); 577 PetscStackCallBLAS("LAPACKsytrs",LAPACKsytrs_("L",&m,&one,mat->v,&mat->lda,mat->pivots,y,&m,&info)); 578 ierr = PetscFPTrapPop();CHKERRQ(ierr); 579 if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"SYTRS Bad solve"); 580 } 581 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix must be factored to solve"); 582 ierr = VecRestoreArrayRead(xx,&x);CHKERRQ(ierr); 583 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 584 ierr = PetscLogFlops(2.0*A->cmap->n*A->cmap->n - A->cmap->n);CHKERRQ(ierr); 585 PetscFunctionReturn(0); 586 } 587 588 /* ---------------------------------------------------------------*/ 589 /* COMMENT: I have chosen to hide row permutation in the pivots, 590 rather than put it in the Mat->row slot.*/ 591 PetscErrorCode MatLUFactor_SeqDense(Mat A,IS row,IS col,const MatFactorInfo *minfo) 592 { 593 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 594 PetscErrorCode ierr; 595 PetscBLASInt n,m,info; 596 597 PetscFunctionBegin; 598 ierr = PetscBLASIntCast(A->cmap->n,&n);CHKERRQ(ierr); 599 ierr = PetscBLASIntCast(A->rmap->n,&m);CHKERRQ(ierr); 600 if (!mat->pivots) { 601 ierr = PetscMalloc1(A->rmap->n,&mat->pivots);CHKERRQ(ierr); 602 ierr = PetscLogObjectMemory((PetscObject)A,A->rmap->n*sizeof(PetscBLASInt));CHKERRQ(ierr); 603 } 604 if (!A->rmap->n || !A->cmap->n) PetscFunctionReturn(0); 605 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr); 606 PetscStackCallBLAS("LAPACKgetrf",LAPACKgetrf_(&m,&n,mat->v,&mat->lda,mat->pivots,&info)); 607 ierr = PetscFPTrapPop();CHKERRQ(ierr); 608 609 if (info<0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Bad argument to LU factorization"); 610 if (info>0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_MAT_LU_ZRPVT,"Bad LU factorization"); 611 612 A->ops->solve = MatSolve_SeqDense; 613 A->ops->matsolve = MatMatSolve_SeqDense; 614 A->ops->solvetranspose = MatSolveTranspose_SeqDense; 615 A->factortype = MAT_FACTOR_LU; 616 617 ierr = PetscFree(A->solvertype);CHKERRQ(ierr); 618 ierr = PetscStrallocpy(MATSOLVERPETSC,&A->solvertype);CHKERRQ(ierr); 619 620 ierr = PetscLogFlops((2.0*A->cmap->n*A->cmap->n*A->cmap->n)/3);CHKERRQ(ierr); 621 PetscFunctionReturn(0); 622 } 623 624 /* Cholesky as L*L^T or L*D*L^T and the symmetric/hermitian complex variants */ 625 PetscErrorCode MatCholeskyFactor_SeqDense(Mat A,IS perm,const MatFactorInfo *factinfo) 626 { 627 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 628 PetscErrorCode ierr; 629 PetscBLASInt info,n; 630 631 PetscFunctionBegin; 632 ierr = PetscBLASIntCast(A->cmap->n,&n);CHKERRQ(ierr); 633 if (!A->rmap->n || !A->cmap->n) PetscFunctionReturn(0); 634 if (A->spd) { 635 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr); 636 PetscStackCallBLAS("LAPACKpotrf",LAPACKpotrf_("L",&n,mat->v,&mat->lda,&info)); 637 ierr = PetscFPTrapPop();CHKERRQ(ierr); 638 #if defined(PETSC_USE_COMPLEX) 639 } else if (A->hermitian) { 640 if (!mat->pivots) { 641 ierr = PetscMalloc1(A->rmap->n,&mat->pivots);CHKERRQ(ierr); 642 ierr = PetscLogObjectMemory((PetscObject)A,A->rmap->n*sizeof(PetscBLASInt));CHKERRQ(ierr); 643 } 644 if (!mat->fwork) { 645 PetscScalar dummy; 646 647 mat->lfwork = -1; 648 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr); 649 PetscStackCallBLAS("LAPACKhetrf",LAPACKhetrf_("L",&n,mat->v,&mat->lda,mat->pivots,&dummy,&mat->lfwork,&info)); 650 ierr = PetscFPTrapPop();CHKERRQ(ierr); 651 mat->lfwork = (PetscInt)PetscRealPart(dummy); 652 ierr = PetscMalloc1(mat->lfwork,&mat->fwork);CHKERRQ(ierr); 653 ierr = PetscLogObjectMemory((PetscObject)A,mat->lfwork*sizeof(PetscBLASInt));CHKERRQ(ierr); 654 } 655 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr); 656 PetscStackCallBLAS("LAPACKhetrf",LAPACKhetrf_("L",&n,mat->v,&mat->lda,mat->pivots,mat->fwork,&mat->lfwork,&info)); 657 ierr = PetscFPTrapPop();CHKERRQ(ierr); 658 #endif 659 } else { /* symmetric case */ 660 if (!mat->pivots) { 661 ierr = PetscMalloc1(A->rmap->n,&mat->pivots);CHKERRQ(ierr); 662 ierr = PetscLogObjectMemory((PetscObject)A,A->rmap->n*sizeof(PetscBLASInt));CHKERRQ(ierr); 663 } 664 if (!mat->fwork) { 665 PetscScalar dummy; 666 667 mat->lfwork = -1; 668 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr); 669 PetscStackCallBLAS("LAPACKsytrf",LAPACKsytrf_("L",&n,mat->v,&mat->lda,mat->pivots,&dummy,&mat->lfwork,&info)); 670 ierr = PetscFPTrapPop();CHKERRQ(ierr); 671 mat->lfwork = (PetscInt)PetscRealPart(dummy); 672 ierr = PetscMalloc1(mat->lfwork,&mat->fwork);CHKERRQ(ierr); 673 ierr = PetscLogObjectMemory((PetscObject)A,mat->lfwork*sizeof(PetscBLASInt));CHKERRQ(ierr); 674 } 675 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr); 676 PetscStackCallBLAS("LAPACKsytrf",LAPACKsytrf_("L",&n,mat->v,&mat->lda,mat->pivots,mat->fwork,&mat->lfwork,&info)); 677 ierr = PetscFPTrapPop();CHKERRQ(ierr); 678 } 679 if (info) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_MAT_CH_ZRPVT,"Bad factorization: zero pivot in row %D",(PetscInt)info-1); 680 681 A->ops->solve = MatSolve_SeqDense; 682 A->ops->matsolve = MatMatSolve_SeqDense; 683 A->ops->solvetranspose = MatSolveTranspose_SeqDense; 684 A->factortype = MAT_FACTOR_CHOLESKY; 685 686 ierr = PetscFree(A->solvertype);CHKERRQ(ierr); 687 ierr = PetscStrallocpy(MATSOLVERPETSC,&A->solvertype);CHKERRQ(ierr); 688 689 ierr = PetscLogFlops((1.0*A->cmap->n*A->cmap->n*A->cmap->n)/3.0);CHKERRQ(ierr); 690 PetscFunctionReturn(0); 691 } 692 693 694 PetscErrorCode MatCholeskyFactorNumeric_SeqDense(Mat fact,Mat A,const MatFactorInfo *info_dummy) 695 { 696 PetscErrorCode ierr; 697 MatFactorInfo info; 698 699 PetscFunctionBegin; 700 info.fill = 1.0; 701 702 ierr = MatDuplicateNoCreate_SeqDense(fact,A,MAT_COPY_VALUES);CHKERRQ(ierr); 703 ierr = (*fact->ops->choleskyfactor)(fact,0,&info);CHKERRQ(ierr); 704 PetscFunctionReturn(0); 705 } 706 707 PetscErrorCode MatCholeskyFactorSymbolic_SeqDense(Mat fact,Mat A,IS row,const MatFactorInfo *info) 708 { 709 PetscFunctionBegin; 710 fact->assembled = PETSC_TRUE; 711 fact->preallocated = PETSC_TRUE; 712 fact->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqDense; 713 fact->ops->solve = MatSolve_SeqDense; 714 fact->ops->matsolve = MatMatSolve_SeqDense; 715 fact->ops->solvetranspose = MatSolveTranspose_SeqDense; 716 PetscFunctionReturn(0); 717 } 718 719 PetscErrorCode MatLUFactorSymbolic_SeqDense(Mat fact,Mat A,IS row,IS col,const MatFactorInfo *info) 720 { 721 PetscFunctionBegin; 722 fact->preallocated = PETSC_TRUE; 723 fact->assembled = PETSC_TRUE; 724 fact->ops->lufactornumeric = MatLUFactorNumeric_SeqDense; 725 fact->ops->solve = MatSolve_SeqDense; 726 fact->ops->matsolve = MatMatSolve_SeqDense; 727 fact->ops->solvetranspose = MatSolveTranspose_SeqDense; 728 PetscFunctionReturn(0); 729 } 730 731 /* uses LAPACK */ 732 PETSC_INTERN PetscErrorCode MatGetFactor_seqdense_petsc(Mat A,MatFactorType ftype,Mat *fact) 733 { 734 PetscErrorCode ierr; 735 736 PetscFunctionBegin; 737 ierr = MatCreate(PetscObjectComm((PetscObject)A),fact);CHKERRQ(ierr); 738 ierr = MatSetSizes(*fact,A->rmap->n,A->cmap->n,A->rmap->n,A->cmap->n);CHKERRQ(ierr); 739 ierr = MatSetType(*fact,MATDENSE);CHKERRQ(ierr); 740 if (ftype == MAT_FACTOR_LU) { 741 (*fact)->ops->lufactorsymbolic = MatLUFactorSymbolic_SeqDense; 742 } else { 743 (*fact)->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_SeqDense; 744 } 745 (*fact)->factortype = ftype; 746 747 ierr = PetscFree((*fact)->solvertype);CHKERRQ(ierr); 748 ierr = PetscStrallocpy(MATSOLVERPETSC,&(*fact)->solvertype);CHKERRQ(ierr); 749 PetscFunctionReturn(0); 750 } 751 752 /* ------------------------------------------------------------------*/ 753 static PetscErrorCode MatSOR_SeqDense(Mat A,Vec bb,PetscReal omega,MatSORType flag,PetscReal shift,PetscInt its,PetscInt lits,Vec xx) 754 { 755 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 756 PetscScalar *x,*v = mat->v,zero = 0.0,xt; 757 const PetscScalar *b; 758 PetscErrorCode ierr; 759 PetscInt m = A->rmap->n,i; 760 PetscBLASInt o = 1,bm; 761 762 PetscFunctionBegin; 763 #if defined(PETSC_HAVE_CUDA) 764 if (A->offloadmask == PETSC_OFFLOAD_GPU) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Not implemented"); 765 #endif 766 if (shift == -1) shift = 0.0; /* negative shift indicates do not error on zero diagonal; this code never zeros on zero diagonal */ 767 ierr = PetscBLASIntCast(m,&bm);CHKERRQ(ierr); 768 if (flag & SOR_ZERO_INITIAL_GUESS) { 769 /* this is a hack fix, should have another version without the second BLASdotu */ 770 ierr = VecSet(xx,zero);CHKERRQ(ierr); 771 } 772 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 773 ierr = VecGetArrayRead(bb,&b);CHKERRQ(ierr); 774 its = its*lits; 775 if (its <= 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Relaxation requires global its %D and local its %D both positive",its,lits); 776 while (its--) { 777 if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP) { 778 for (i=0; i<m; i++) { 779 PetscStackCallBLAS("BLASdotu",xt = b[i] - BLASdotu_(&bm,v+i,&bm,x,&o)); 780 x[i] = (1. - omega)*x[i] + omega*(xt+v[i + i*m]*x[i])/(v[i + i*m]+shift); 781 } 782 } 783 if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP) { 784 for (i=m-1; i>=0; i--) { 785 PetscStackCallBLAS("BLASdotu",xt = b[i] - BLASdotu_(&bm,v+i,&bm,x,&o)); 786 x[i] = (1. - omega)*x[i] + omega*(xt+v[i + i*m]*x[i])/(v[i + i*m]+shift); 787 } 788 } 789 } 790 ierr = VecRestoreArrayRead(bb,&b);CHKERRQ(ierr); 791 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 792 PetscFunctionReturn(0); 793 } 794 795 /* -----------------------------------------------------------------*/ 796 PetscErrorCode MatMultTranspose_SeqDense(Mat A,Vec xx,Vec yy) 797 { 798 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 799 const PetscScalar *v = mat->v,*x; 800 PetscScalar *y; 801 PetscErrorCode ierr; 802 PetscBLASInt m, n,_One=1; 803 PetscScalar _DOne=1.0,_DZero=0.0; 804 805 PetscFunctionBegin; 806 ierr = PetscBLASIntCast(A->rmap->n,&m);CHKERRQ(ierr); 807 ierr = PetscBLASIntCast(A->cmap->n,&n);CHKERRQ(ierr); 808 ierr = VecGetArrayRead(xx,&x);CHKERRQ(ierr); 809 ierr = VecGetArrayWrite(yy,&y);CHKERRQ(ierr); 810 if (!A->rmap->n || !A->cmap->n) { 811 PetscBLASInt i; 812 for (i=0; i<n; i++) y[i] = 0.0; 813 } else { 814 PetscStackCallBLAS("BLASgemv",BLASgemv_("T",&m,&n,&_DOne,v,&mat->lda,x,&_One,&_DZero,y,&_One)); 815 ierr = PetscLogFlops(2.0*A->rmap->n*A->cmap->n - A->cmap->n);CHKERRQ(ierr); 816 } 817 ierr = VecRestoreArrayRead(xx,&x);CHKERRQ(ierr); 818 ierr = VecRestoreArrayWrite(yy,&y);CHKERRQ(ierr); 819 PetscFunctionReturn(0); 820 } 821 822 PetscErrorCode MatMult_SeqDense(Mat A,Vec xx,Vec yy) 823 { 824 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 825 PetscScalar *y,_DOne=1.0,_DZero=0.0; 826 PetscErrorCode ierr; 827 PetscBLASInt m, n, _One=1; 828 const PetscScalar *v = mat->v,*x; 829 830 PetscFunctionBegin; 831 ierr = PetscBLASIntCast(A->rmap->n,&m);CHKERRQ(ierr); 832 ierr = PetscBLASIntCast(A->cmap->n,&n);CHKERRQ(ierr); 833 ierr = VecGetArrayRead(xx,&x);CHKERRQ(ierr); 834 ierr = VecGetArrayWrite(yy,&y);CHKERRQ(ierr); 835 if (!A->rmap->n || !A->cmap->n) { 836 PetscBLASInt i; 837 for (i=0; i<m; i++) y[i] = 0.0; 838 } else { 839 PetscStackCallBLAS("BLASgemv",BLASgemv_("N",&m,&n,&_DOne,v,&(mat->lda),x,&_One,&_DZero,y,&_One)); 840 ierr = PetscLogFlops(2.0*A->rmap->n*A->cmap->n - A->rmap->n);CHKERRQ(ierr); 841 } 842 ierr = VecRestoreArrayRead(xx,&x);CHKERRQ(ierr); 843 ierr = VecRestoreArrayWrite(yy,&y);CHKERRQ(ierr); 844 PetscFunctionReturn(0); 845 } 846 847 PetscErrorCode MatMultAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy) 848 { 849 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 850 const PetscScalar *v = mat->v,*x; 851 PetscScalar *y,_DOne=1.0; 852 PetscErrorCode ierr; 853 PetscBLASInt m, n, _One=1; 854 855 PetscFunctionBegin; 856 ierr = PetscBLASIntCast(A->rmap->n,&m);CHKERRQ(ierr); 857 ierr = PetscBLASIntCast(A->cmap->n,&n);CHKERRQ(ierr); 858 if (!A->rmap->n || !A->cmap->n) PetscFunctionReturn(0); 859 if (zz != yy) {ierr = VecCopy(zz,yy);CHKERRQ(ierr);} 860 ierr = VecGetArrayRead(xx,&x);CHKERRQ(ierr); 861 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 862 PetscStackCallBLAS("BLASgemv",BLASgemv_("N",&m,&n,&_DOne,v,&(mat->lda),x,&_One,&_DOne,y,&_One)); 863 ierr = VecRestoreArrayRead(xx,&x);CHKERRQ(ierr); 864 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 865 ierr = PetscLogFlops(2.0*A->rmap->n*A->cmap->n);CHKERRQ(ierr); 866 PetscFunctionReturn(0); 867 } 868 869 PetscErrorCode MatMultTransposeAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy) 870 { 871 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 872 const PetscScalar *v = mat->v,*x; 873 PetscScalar *y; 874 PetscErrorCode ierr; 875 PetscBLASInt m, n, _One=1; 876 PetscScalar _DOne=1.0; 877 878 PetscFunctionBegin; 879 ierr = PetscBLASIntCast(A->rmap->n,&m);CHKERRQ(ierr); 880 ierr = PetscBLASIntCast(A->cmap->n,&n);CHKERRQ(ierr); 881 if (!A->rmap->n || !A->cmap->n) PetscFunctionReturn(0); 882 if (zz != yy) {ierr = VecCopy(zz,yy);CHKERRQ(ierr);} 883 ierr = VecGetArrayRead(xx,&x);CHKERRQ(ierr); 884 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 885 PetscStackCallBLAS("BLASgemv",BLASgemv_("T",&m,&n,&_DOne,v,&(mat->lda),x,&_One,&_DOne,y,&_One)); 886 ierr = VecRestoreArrayRead(xx,&x);CHKERRQ(ierr); 887 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 888 ierr = PetscLogFlops(2.0*A->rmap->n*A->cmap->n);CHKERRQ(ierr); 889 PetscFunctionReturn(0); 890 } 891 892 /* -----------------------------------------------------------------*/ 893 static PetscErrorCode MatGetRow_SeqDense(Mat A,PetscInt row,PetscInt *ncols,PetscInt **cols,PetscScalar **vals) 894 { 895 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 896 PetscErrorCode ierr; 897 PetscInt i; 898 899 PetscFunctionBegin; 900 *ncols = A->cmap->n; 901 if (cols) { 902 ierr = PetscMalloc1(A->cmap->n+1,cols);CHKERRQ(ierr); 903 for (i=0; i<A->cmap->n; i++) (*cols)[i] = i; 904 } 905 if (vals) { 906 const PetscScalar *v; 907 908 ierr = MatDenseGetArrayRead(A,&v);CHKERRQ(ierr); 909 ierr = PetscMalloc1(A->cmap->n+1,vals);CHKERRQ(ierr); 910 v += row; 911 for (i=0; i<A->cmap->n; i++) {(*vals)[i] = *v; v += mat->lda;} 912 ierr = MatDenseRestoreArrayRead(A,&v);CHKERRQ(ierr); 913 } 914 PetscFunctionReturn(0); 915 } 916 917 static PetscErrorCode MatRestoreRow_SeqDense(Mat A,PetscInt row,PetscInt *ncols,PetscInt **cols,PetscScalar **vals) 918 { 919 PetscErrorCode ierr; 920 921 PetscFunctionBegin; 922 if (cols) {ierr = PetscFree(*cols);CHKERRQ(ierr);} 923 if (vals) {ierr = PetscFree(*vals);CHKERRQ(ierr); } 924 PetscFunctionReturn(0); 925 } 926 /* ----------------------------------------------------------------*/ 927 static PetscErrorCode MatSetValues_SeqDense(Mat A,PetscInt m,const PetscInt indexm[],PetscInt n,const PetscInt indexn[],const PetscScalar v[],InsertMode addv) 928 { 929 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 930 PetscScalar *av; 931 PetscInt i,j,idx=0; 932 #if defined(PETSC_HAVE_CUDA) 933 PetscOffloadMask oldf; 934 #endif 935 PetscErrorCode ierr; 936 937 PetscFunctionBegin; 938 ierr = MatDenseGetArray(A,&av);CHKERRQ(ierr); 939 if (!mat->roworiented) { 940 if (addv == INSERT_VALUES) { 941 for (j=0; j<n; j++) { 942 if (indexn[j] < 0) {idx += m; continue;} 943 if (PetscUnlikelyDebug(indexn[j] >= A->cmap->n)) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %D max %D",indexn[j],A->cmap->n-1); 944 for (i=0; i<m; i++) { 945 if (indexm[i] < 0) {idx++; continue;} 946 if (PetscUnlikelyDebug(indexm[i] >= A->rmap->n)) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",indexm[i],A->rmap->n-1); 947 av[indexn[j]*mat->lda + indexm[i]] = v[idx++]; 948 } 949 } 950 } else { 951 for (j=0; j<n; j++) { 952 if (indexn[j] < 0) {idx += m; continue;} 953 if (PetscUnlikelyDebug(indexn[j] >= A->cmap->n)) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %D max %D",indexn[j],A->cmap->n-1); 954 for (i=0; i<m; i++) { 955 if (indexm[i] < 0) {idx++; continue;} 956 if (PetscUnlikelyDebug(indexm[i] >= A->rmap->n)) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",indexm[i],A->rmap->n-1); 957 av[indexn[j]*mat->lda + indexm[i]] += v[idx++]; 958 } 959 } 960 } 961 } else { 962 if (addv == INSERT_VALUES) { 963 for (i=0; i<m; i++) { 964 if (indexm[i] < 0) { idx += n; continue;} 965 if (PetscUnlikelyDebug(indexm[i] >= A->rmap->n)) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",indexm[i],A->rmap->n-1); 966 for (j=0; j<n; j++) { 967 if (indexn[j] < 0) { idx++; continue;} 968 if (PetscUnlikelyDebug(indexn[j] >= A->cmap->n)) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %D max %D",indexn[j],A->cmap->n-1); 969 av[indexn[j]*mat->lda + indexm[i]] = v[idx++]; 970 } 971 } 972 } else { 973 for (i=0; i<m; i++) { 974 if (indexm[i] < 0) { idx += n; continue;} 975 if (PetscUnlikelyDebug(indexm[i] >= A->rmap->n)) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",indexm[i],A->rmap->n-1); 976 for (j=0; j<n; j++) { 977 if (indexn[j] < 0) { idx++; continue;} 978 if (PetscUnlikelyDebug(indexn[j] >= A->cmap->n)) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %D max %D",indexn[j],A->cmap->n-1); 979 av[indexn[j]*mat->lda + indexm[i]] += v[idx++]; 980 } 981 } 982 } 983 } 984 /* hack to prevent unneeded copy to the GPU while returning the array */ 985 #if defined(PETSC_HAVE_CUDA) 986 oldf = A->offloadmask; 987 A->offloadmask = PETSC_OFFLOAD_GPU; 988 #endif 989 ierr = MatDenseRestoreArray(A,&av);CHKERRQ(ierr); 990 #if defined(PETSC_HAVE_CUDA) 991 A->offloadmask = (oldf == PETSC_OFFLOAD_UNALLOCATED ? PETSC_OFFLOAD_UNALLOCATED : PETSC_OFFLOAD_CPU); 992 #endif 993 PetscFunctionReturn(0); 994 } 995 996 static PetscErrorCode MatGetValues_SeqDense(Mat A,PetscInt m,const PetscInt indexm[],PetscInt n,const PetscInt indexn[],PetscScalar v[]) 997 { 998 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 999 const PetscScalar *vv; 1000 PetscInt i,j; 1001 PetscErrorCode ierr; 1002 1003 PetscFunctionBegin; 1004 ierr = MatDenseGetArrayRead(A,&vv);CHKERRQ(ierr); 1005 /* row-oriented output */ 1006 for (i=0; i<m; i++) { 1007 if (indexm[i] < 0) {v += n;continue;} 1008 if (indexm[i] >= A->rmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row %D requested larger than number rows %D",indexm[i],A->rmap->n); 1009 for (j=0; j<n; j++) { 1010 if (indexn[j] < 0) {v++; continue;} 1011 if (indexn[j] >= A->cmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Column %D requested larger than number columns %D",indexn[j],A->cmap->n); 1012 *v++ = vv[indexn[j]*mat->lda + indexm[i]]; 1013 } 1014 } 1015 ierr = MatDenseRestoreArrayRead(A,&vv);CHKERRQ(ierr); 1016 PetscFunctionReturn(0); 1017 } 1018 1019 /* -----------------------------------------------------------------*/ 1020 1021 PetscErrorCode MatView_Dense_Binary(Mat mat,PetscViewer viewer) 1022 { 1023 PetscErrorCode ierr; 1024 PetscBool skipHeader; 1025 PetscViewerFormat format; 1026 PetscInt header[4],M,N,m,lda,i,j,k; 1027 const PetscScalar *v; 1028 PetscScalar *vwork; 1029 1030 PetscFunctionBegin; 1031 ierr = PetscViewerSetUp(viewer);CHKERRQ(ierr); 1032 ierr = PetscViewerBinaryGetSkipHeader(viewer,&skipHeader);CHKERRQ(ierr); 1033 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 1034 if (skipHeader) format = PETSC_VIEWER_NATIVE; 1035 1036 ierr = MatGetSize(mat,&M,&N);CHKERRQ(ierr); 1037 1038 /* write matrix header */ 1039 header[0] = MAT_FILE_CLASSID; header[1] = M; header[2] = N; 1040 header[3] = (format == PETSC_VIEWER_NATIVE) ? MATRIX_BINARY_FORMAT_DENSE : M*N; 1041 if (!skipHeader) {ierr = PetscViewerBinaryWrite(viewer,header,4,PETSC_INT);CHKERRQ(ierr);} 1042 1043 ierr = MatGetLocalSize(mat,&m,NULL);CHKERRQ(ierr); 1044 if (format != PETSC_VIEWER_NATIVE) { 1045 PetscInt nnz = m*N, *iwork; 1046 /* store row lengths for each row */ 1047 ierr = PetscMalloc1(nnz,&iwork);CHKERRQ(ierr); 1048 for (i=0; i<m; i++) iwork[i] = N; 1049 ierr = PetscViewerBinaryWriteAll(viewer,iwork,m,PETSC_DETERMINE,PETSC_DETERMINE,PETSC_INT);CHKERRQ(ierr); 1050 /* store column indices (zero start index) */ 1051 for (k=0, i=0; i<m; i++) 1052 for (j=0; j<N; j++, k++) 1053 iwork[k] = j; 1054 ierr = PetscViewerBinaryWriteAll(viewer,iwork,nnz,PETSC_DETERMINE,PETSC_DETERMINE,PETSC_INT);CHKERRQ(ierr); 1055 ierr = PetscFree(iwork);CHKERRQ(ierr); 1056 } 1057 /* store matrix values as a dense matrix in row major order */ 1058 ierr = PetscMalloc1(m*N,&vwork);CHKERRQ(ierr); 1059 ierr = MatDenseGetArrayRead(mat,&v);CHKERRQ(ierr); 1060 ierr = MatDenseGetLDA(mat,&lda);CHKERRQ(ierr); 1061 for (k=0, i=0; i<m; i++) 1062 for (j=0; j<N; j++, k++) 1063 vwork[k] = v[i+lda*j]; 1064 ierr = MatDenseRestoreArrayRead(mat,&v);CHKERRQ(ierr); 1065 ierr = PetscViewerBinaryWriteAll(viewer,vwork,m*N,PETSC_DETERMINE,PETSC_DETERMINE,PETSC_SCALAR);CHKERRQ(ierr); 1066 ierr = PetscFree(vwork);CHKERRQ(ierr); 1067 PetscFunctionReturn(0); 1068 } 1069 1070 PetscErrorCode MatLoad_Dense_Binary(Mat mat,PetscViewer viewer) 1071 { 1072 PetscErrorCode ierr; 1073 PetscBool skipHeader; 1074 PetscInt header[4],M,N,m,nz,lda,i,j,k; 1075 PetscInt rows,cols; 1076 PetscScalar *v,*vwork; 1077 1078 PetscFunctionBegin; 1079 ierr = PetscViewerSetUp(viewer);CHKERRQ(ierr); 1080 ierr = PetscViewerBinaryGetSkipHeader(viewer,&skipHeader);CHKERRQ(ierr); 1081 1082 if (!skipHeader) { 1083 ierr = PetscViewerBinaryRead(viewer,header,4,NULL,PETSC_INT);CHKERRQ(ierr); 1084 if (header[0] != MAT_FILE_CLASSID) SETERRQ(PetscObjectComm((PetscObject)viewer),PETSC_ERR_FILE_UNEXPECTED,"Not a matrix object in file"); 1085 M = header[1]; N = header[2]; 1086 if (M < 0) SETERRQ1(PetscObjectComm((PetscObject)viewer),PETSC_ERR_FILE_UNEXPECTED,"Matrix row size (%D) in file is negative",M); 1087 if (N < 0) SETERRQ1(PetscObjectComm((PetscObject)viewer),PETSC_ERR_FILE_UNEXPECTED,"Matrix column size (%D) in file is negative",N); 1088 nz = header[3]; 1089 if (nz != MATRIX_BINARY_FORMAT_DENSE && nz < 0) SETERRQ1(PetscObjectComm((PetscObject)viewer),PETSC_ERR_FILE_UNEXPECTED,"Unknown matrix format %D in file",nz); 1090 } else { 1091 ierr = MatGetSize(mat,&M,&N);CHKERRQ(ierr); 1092 if (M < 0 || N < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Matrix binary file header was skipped, thus the user must specify the global sizes of input matrix"); 1093 nz = MATRIX_BINARY_FORMAT_DENSE; 1094 } 1095 1096 /* setup global sizes if not set */ 1097 if (mat->rmap->N < 0) mat->rmap->N = M; 1098 if (mat->cmap->N < 0) mat->cmap->N = N; 1099 ierr = MatSetUp(mat);CHKERRQ(ierr); 1100 /* check if global sizes are correct */ 1101 ierr = MatGetSize(mat,&rows,&cols);CHKERRQ(ierr); 1102 if (M != rows || N != cols) SETERRQ4(PetscObjectComm((PetscObject)viewer),PETSC_ERR_FILE_UNEXPECTED, "Matrix in file of different sizes (%d, %d) than the input matrix (%d, %d)",M,N,rows,cols); 1103 1104 ierr = MatGetSize(mat,NULL,&N);CHKERRQ(ierr); 1105 ierr = MatGetLocalSize(mat,&m,NULL);CHKERRQ(ierr); 1106 ierr = MatDenseGetArray(mat,&v);CHKERRQ(ierr); 1107 ierr = MatDenseGetLDA(mat,&lda);CHKERRQ(ierr); 1108 if (nz == MATRIX_BINARY_FORMAT_DENSE) { /* matrix in file is dense format */ 1109 PetscInt nnz = m*N; 1110 /* read in matrix values */ 1111 ierr = PetscMalloc1(nnz,&vwork);CHKERRQ(ierr); 1112 ierr = PetscViewerBinaryReadAll(viewer,vwork,nnz,PETSC_DETERMINE,PETSC_DETERMINE,PETSC_SCALAR);CHKERRQ(ierr); 1113 /* store values in column major order */ 1114 for (j=0; j<N; j++) 1115 for (i=0; i<m; i++) 1116 v[i+lda*j] = vwork[i*N+j]; 1117 ierr = PetscFree(vwork);CHKERRQ(ierr); 1118 } else { /* matrix in file is sparse format */ 1119 PetscInt nnz = 0, *rlens, *icols; 1120 /* read in row lengths */ 1121 ierr = PetscMalloc1(m,&rlens);CHKERRQ(ierr); 1122 ierr = PetscViewerBinaryReadAll(viewer,rlens,m,PETSC_DETERMINE,PETSC_DETERMINE,PETSC_INT);CHKERRQ(ierr); 1123 for (i=0; i<m; i++) nnz += rlens[i]; 1124 /* read in column indices and values */ 1125 ierr = PetscMalloc2(nnz,&icols,nnz,&vwork);CHKERRQ(ierr); 1126 ierr = PetscViewerBinaryReadAll(viewer,icols,nnz,PETSC_DETERMINE,PETSC_DETERMINE,PETSC_INT);CHKERRQ(ierr); 1127 ierr = PetscViewerBinaryReadAll(viewer,vwork,nnz,PETSC_DETERMINE,PETSC_DETERMINE,PETSC_SCALAR);CHKERRQ(ierr); 1128 /* store values in column major order */ 1129 for (k=0, i=0; i<m; i++) 1130 for (j=0; j<rlens[i]; j++, k++) 1131 v[i+lda*icols[k]] = vwork[k]; 1132 ierr = PetscFree(rlens);CHKERRQ(ierr); 1133 ierr = PetscFree2(icols,vwork);CHKERRQ(ierr); 1134 } 1135 ierr = MatDenseRestoreArray(mat,&v);CHKERRQ(ierr); 1136 ierr = MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1137 ierr = MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1138 PetscFunctionReturn(0); 1139 } 1140 1141 PetscErrorCode MatLoad_SeqDense(Mat newMat, PetscViewer viewer) 1142 { 1143 PetscBool isbinary, ishdf5; 1144 PetscErrorCode ierr; 1145 1146 PetscFunctionBegin; 1147 PetscValidHeaderSpecific(newMat,MAT_CLASSID,1); 1148 PetscValidHeaderSpecific(viewer,PETSC_VIEWER_CLASSID,2); 1149 /* force binary viewer to load .info file if it has not yet done so */ 1150 ierr = PetscViewerSetUp(viewer);CHKERRQ(ierr); 1151 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);CHKERRQ(ierr); 1152 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERHDF5, &ishdf5);CHKERRQ(ierr); 1153 if (isbinary) { 1154 ierr = MatLoad_Dense_Binary(newMat,viewer);CHKERRQ(ierr); 1155 } else if (ishdf5) { 1156 #if defined(PETSC_HAVE_HDF5) 1157 ierr = MatLoad_Dense_HDF5(newMat,viewer);CHKERRQ(ierr); 1158 #else 1159 SETERRQ(PetscObjectComm((PetscObject)newMat),PETSC_ERR_SUP,"HDF5 not supported in this build.\nPlease reconfigure using --download-hdf5"); 1160 #endif 1161 } else { 1162 SETERRQ2(PetscObjectComm((PetscObject)newMat),PETSC_ERR_SUP,"Viewer type %s not yet supported for reading %s matrices",((PetscObject)viewer)->type_name,((PetscObject)newMat)->type_name); 1163 } 1164 PetscFunctionReturn(0); 1165 } 1166 1167 static PetscErrorCode MatView_SeqDense_ASCII(Mat A,PetscViewer viewer) 1168 { 1169 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1170 PetscErrorCode ierr; 1171 PetscInt i,j; 1172 const char *name; 1173 PetscScalar *v,*av; 1174 PetscViewerFormat format; 1175 #if defined(PETSC_USE_COMPLEX) 1176 PetscBool allreal = PETSC_TRUE; 1177 #endif 1178 1179 PetscFunctionBegin; 1180 ierr = MatDenseGetArrayRead(A,(const PetscScalar**)&av);CHKERRQ(ierr); 1181 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 1182 if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) { 1183 PetscFunctionReturn(0); /* do nothing for now */ 1184 } else if (format == PETSC_VIEWER_ASCII_COMMON) { 1185 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);CHKERRQ(ierr); 1186 for (i=0; i<A->rmap->n; i++) { 1187 v = av + i; 1188 ierr = PetscViewerASCIIPrintf(viewer,"row %D:",i);CHKERRQ(ierr); 1189 for (j=0; j<A->cmap->n; j++) { 1190 #if defined(PETSC_USE_COMPLEX) 1191 if (PetscRealPart(*v) != 0.0 && PetscImaginaryPart(*v) != 0.0) { 1192 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g + %g i) ",j,(double)PetscRealPart(*v),(double)PetscImaginaryPart(*v));CHKERRQ(ierr); 1193 } else if (PetscRealPart(*v)) { 1194 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g) ",j,(double)PetscRealPart(*v));CHKERRQ(ierr); 1195 } 1196 #else 1197 if (*v) { 1198 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g) ",j,(double)*v);CHKERRQ(ierr); 1199 } 1200 #endif 1201 v += a->lda; 1202 } 1203 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 1204 } 1205 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);CHKERRQ(ierr); 1206 } else { 1207 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);CHKERRQ(ierr); 1208 #if defined(PETSC_USE_COMPLEX) 1209 /* determine if matrix has all real values */ 1210 v = av; 1211 for (i=0; i<A->rmap->n*A->cmap->n; i++) { 1212 if (PetscImaginaryPart(v[i])) { allreal = PETSC_FALSE; break;} 1213 } 1214 #endif 1215 if (format == PETSC_VIEWER_ASCII_MATLAB) { 1216 ierr = PetscObjectGetName((PetscObject)A,&name);CHKERRQ(ierr); 1217 ierr = PetscViewerASCIIPrintf(viewer,"%% Size = %D %D \n",A->rmap->n,A->cmap->n);CHKERRQ(ierr); 1218 ierr = PetscViewerASCIIPrintf(viewer,"%s = zeros(%D,%D);\n",name,A->rmap->n,A->cmap->n);CHKERRQ(ierr); 1219 ierr = PetscViewerASCIIPrintf(viewer,"%s = [\n",name);CHKERRQ(ierr); 1220 } 1221 1222 for (i=0; i<A->rmap->n; i++) { 1223 v = av + i; 1224 for (j=0; j<A->cmap->n; j++) { 1225 #if defined(PETSC_USE_COMPLEX) 1226 if (allreal) { 1227 ierr = PetscViewerASCIIPrintf(viewer,"%18.16e ",(double)PetscRealPart(*v));CHKERRQ(ierr); 1228 } else { 1229 ierr = PetscViewerASCIIPrintf(viewer,"%18.16e + %18.16ei ",(double)PetscRealPart(*v),(double)PetscImaginaryPart(*v));CHKERRQ(ierr); 1230 } 1231 #else 1232 ierr = PetscViewerASCIIPrintf(viewer,"%18.16e ",(double)*v);CHKERRQ(ierr); 1233 #endif 1234 v += a->lda; 1235 } 1236 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 1237 } 1238 if (format == PETSC_VIEWER_ASCII_MATLAB) { 1239 ierr = PetscViewerASCIIPrintf(viewer,"];\n");CHKERRQ(ierr); 1240 } 1241 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);CHKERRQ(ierr); 1242 } 1243 ierr = MatDenseRestoreArrayRead(A,(const PetscScalar**)&av);CHKERRQ(ierr); 1244 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 1245 PetscFunctionReturn(0); 1246 } 1247 1248 #include <petscdraw.h> 1249 static PetscErrorCode MatView_SeqDense_Draw_Zoom(PetscDraw draw,void *Aa) 1250 { 1251 Mat A = (Mat) Aa; 1252 PetscErrorCode ierr; 1253 PetscInt m = A->rmap->n,n = A->cmap->n,i,j; 1254 int color = PETSC_DRAW_WHITE; 1255 const PetscScalar *v; 1256 PetscViewer viewer; 1257 PetscReal xl,yl,xr,yr,x_l,x_r,y_l,y_r; 1258 PetscViewerFormat format; 1259 1260 PetscFunctionBegin; 1261 ierr = PetscObjectQuery((PetscObject)A,"Zoomviewer",(PetscObject*)&viewer);CHKERRQ(ierr); 1262 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 1263 ierr = PetscDrawGetCoordinates(draw,&xl,&yl,&xr,&yr);CHKERRQ(ierr); 1264 1265 /* Loop over matrix elements drawing boxes */ 1266 ierr = MatDenseGetArrayRead(A,&v);CHKERRQ(ierr); 1267 if (format != PETSC_VIEWER_DRAW_CONTOUR) { 1268 ierr = PetscDrawCollectiveBegin(draw);CHKERRQ(ierr); 1269 /* Blue for negative and Red for positive */ 1270 for (j = 0; j < n; j++) { 1271 x_l = j; x_r = x_l + 1.0; 1272 for (i = 0; i < m; i++) { 1273 y_l = m - i - 1.0; 1274 y_r = y_l + 1.0; 1275 if (PetscRealPart(v[j*m+i]) > 0.) color = PETSC_DRAW_RED; 1276 else if (PetscRealPart(v[j*m+i]) < 0.) color = PETSC_DRAW_BLUE; 1277 else continue; 1278 ierr = PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);CHKERRQ(ierr); 1279 } 1280 } 1281 ierr = PetscDrawCollectiveEnd(draw);CHKERRQ(ierr); 1282 } else { 1283 /* use contour shading to indicate magnitude of values */ 1284 /* first determine max of all nonzero values */ 1285 PetscReal minv = 0.0, maxv = 0.0; 1286 PetscDraw popup; 1287 1288 for (i=0; i < m*n; i++) { 1289 if (PetscAbsScalar(v[i]) > maxv) maxv = PetscAbsScalar(v[i]); 1290 } 1291 if (minv >= maxv) maxv = minv + PETSC_SMALL; 1292 ierr = PetscDrawGetPopup(draw,&popup);CHKERRQ(ierr); 1293 ierr = PetscDrawScalePopup(popup,minv,maxv);CHKERRQ(ierr); 1294 1295 ierr = PetscDrawCollectiveBegin(draw);CHKERRQ(ierr); 1296 for (j=0; j<n; j++) { 1297 x_l = j; 1298 x_r = x_l + 1.0; 1299 for (i=0; i<m; i++) { 1300 y_l = m - i - 1.0; 1301 y_r = y_l + 1.0; 1302 color = PetscDrawRealToColor(PetscAbsScalar(v[j*m+i]),minv,maxv); 1303 ierr = PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);CHKERRQ(ierr); 1304 } 1305 } 1306 ierr = PetscDrawCollectiveEnd(draw);CHKERRQ(ierr); 1307 } 1308 ierr = MatDenseRestoreArrayRead(A,&v);CHKERRQ(ierr); 1309 PetscFunctionReturn(0); 1310 } 1311 1312 static PetscErrorCode MatView_SeqDense_Draw(Mat A,PetscViewer viewer) 1313 { 1314 PetscDraw draw; 1315 PetscBool isnull; 1316 PetscReal xr,yr,xl,yl,h,w; 1317 PetscErrorCode ierr; 1318 1319 PetscFunctionBegin; 1320 ierr = PetscViewerDrawGetDraw(viewer,0,&draw);CHKERRQ(ierr); 1321 ierr = PetscDrawIsNull(draw,&isnull);CHKERRQ(ierr); 1322 if (isnull) PetscFunctionReturn(0); 1323 1324 xr = A->cmap->n; yr = A->rmap->n; h = yr/10.0; w = xr/10.0; 1325 xr += w; yr += h; xl = -w; yl = -h; 1326 ierr = PetscDrawSetCoordinates(draw,xl,yl,xr,yr);CHKERRQ(ierr); 1327 ierr = PetscObjectCompose((PetscObject)A,"Zoomviewer",(PetscObject)viewer);CHKERRQ(ierr); 1328 ierr = PetscDrawZoom(draw,MatView_SeqDense_Draw_Zoom,A);CHKERRQ(ierr); 1329 ierr = PetscObjectCompose((PetscObject)A,"Zoomviewer",NULL);CHKERRQ(ierr); 1330 ierr = PetscDrawSave(draw);CHKERRQ(ierr); 1331 PetscFunctionReturn(0); 1332 } 1333 1334 PetscErrorCode MatView_SeqDense(Mat A,PetscViewer viewer) 1335 { 1336 PetscErrorCode ierr; 1337 PetscBool iascii,isbinary,isdraw; 1338 1339 PetscFunctionBegin; 1340 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 1341 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);CHKERRQ(ierr); 1342 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);CHKERRQ(ierr); 1343 1344 if (iascii) { 1345 ierr = MatView_SeqDense_ASCII(A,viewer);CHKERRQ(ierr); 1346 } else if (isbinary) { 1347 ierr = MatView_Dense_Binary(A,viewer);CHKERRQ(ierr); 1348 } else if (isdraw) { 1349 ierr = MatView_SeqDense_Draw(A,viewer);CHKERRQ(ierr); 1350 } 1351 PetscFunctionReturn(0); 1352 } 1353 1354 static PetscErrorCode MatDensePlaceArray_SeqDense(Mat A,const PetscScalar *array) 1355 { 1356 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1357 1358 PetscFunctionBegin; 1359 a->unplacedarray = a->v; 1360 a->unplaced_user_alloc = a->user_alloc; 1361 a->v = (PetscScalar*) array; 1362 a->user_alloc = PETSC_TRUE; 1363 #if defined(PETSC_HAVE_CUDA) 1364 A->offloadmask = PETSC_OFFLOAD_CPU; 1365 #endif 1366 PetscFunctionReturn(0); 1367 } 1368 1369 static PetscErrorCode MatDenseResetArray_SeqDense(Mat A) 1370 { 1371 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1372 1373 PetscFunctionBegin; 1374 a->v = a->unplacedarray; 1375 a->user_alloc = a->unplaced_user_alloc; 1376 a->unplacedarray = NULL; 1377 #if defined(PETSC_HAVE_CUDA) 1378 A->offloadmask = PETSC_OFFLOAD_CPU; 1379 #endif 1380 PetscFunctionReturn(0); 1381 } 1382 1383 PetscErrorCode MatDestroy_SeqDense(Mat mat) 1384 { 1385 Mat_SeqDense *l = (Mat_SeqDense*)mat->data; 1386 PetscErrorCode ierr; 1387 1388 PetscFunctionBegin; 1389 #if defined(PETSC_USE_LOG) 1390 PetscLogObjectState((PetscObject)mat,"Rows %D Cols %D",mat->rmap->n,mat->cmap->n); 1391 #endif 1392 ierr = PetscFree(l->pivots);CHKERRQ(ierr); 1393 ierr = PetscFree(l->fwork);CHKERRQ(ierr); 1394 ierr = MatDestroy(&l->ptapwork);CHKERRQ(ierr); 1395 if (!l->user_alloc) {ierr = PetscFree(l->v);CHKERRQ(ierr);} 1396 if (!l->unplaced_user_alloc) {ierr = PetscFree(l->unplacedarray);CHKERRQ(ierr);} 1397 ierr = PetscFree(mat->data);CHKERRQ(ierr); 1398 1399 ierr = PetscObjectChangeTypeName((PetscObject)mat,0);CHKERRQ(ierr); 1400 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetLDA_C",NULL);CHKERRQ(ierr); 1401 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetArray_C",NULL);CHKERRQ(ierr); 1402 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatDenseRestoreArray_C",NULL);CHKERRQ(ierr); 1403 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatDensePlaceArray_C",NULL);CHKERRQ(ierr); 1404 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatDenseResetArray_C",NULL);CHKERRQ(ierr); 1405 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetArrayRead_C",NULL);CHKERRQ(ierr); 1406 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatDenseRestoreArrayRead_C",NULL);CHKERRQ(ierr); 1407 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatConvert_seqdense_seqaij_C",NULL);CHKERRQ(ierr); 1408 #if defined(PETSC_HAVE_ELEMENTAL) 1409 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatConvert_seqdense_elemental_C",NULL);CHKERRQ(ierr); 1410 #endif 1411 #if defined(PETSC_HAVE_CUDA) 1412 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatConvert_seqdense_seqdensecuda_C",NULL);CHKERRQ(ierr); 1413 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatProductSetFromOptions_seqdensecuda_seqdensecuda_C",NULL);CHKERRQ(ierr); 1414 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatProductSetFromOptions_seqdensecuda_seqdense_C",NULL);CHKERRQ(ierr); 1415 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatProductSetFromOptions_seqaijcusparse_seqdense_C",NULL);CHKERRQ(ierr); 1416 #endif 1417 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatSeqDenseSetPreallocation_C",NULL);CHKERRQ(ierr); 1418 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatProductSetFromOptions_seqaij_seqdense_C",NULL);CHKERRQ(ierr); 1419 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatProductSetFromOptions_seqdense_seqdense_C",NULL);CHKERRQ(ierr); 1420 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatMatMultSymbolic_seqaij_seqdense_C",NULL);CHKERRQ(ierr); 1421 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatMatMultNumeric_seqaij_seqdense_C",NULL);CHKERRQ(ierr); 1422 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatProductSetFromOptions_seqbaij_seqdense_C",NULL);CHKERRQ(ierr); 1423 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatMatMultSymbolic_seqbaij_seqdense_C",NULL);CHKERRQ(ierr); 1424 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatMatMultNumeric_seqbaij_seqdense_C",NULL);CHKERRQ(ierr); 1425 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatProductSetFromOptions_seqsbaij_seqdense_C",NULL);CHKERRQ(ierr); 1426 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatMatMultSymbolic_seqsbaij_seqdense_C",NULL);CHKERRQ(ierr); 1427 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatMatMultNumeric_seqsbaij_seqdense_C",NULL);CHKERRQ(ierr); 1428 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatMatMultSymbolic_seqaijperm_seqdense_C",NULL);CHKERRQ(ierr); 1429 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatMatMultNumeric_seqaijperm_seqdense_C",NULL);CHKERRQ(ierr); 1430 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatMatMultSymbolic_seqaijsell_seqdense_C",NULL);CHKERRQ(ierr); 1431 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatMatMultNumeric_seqaijsell_seqdense_C",NULL);CHKERRQ(ierr); 1432 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatMatMultSymbolic_seqaijmkl_seqdense_C",NULL);CHKERRQ(ierr); 1433 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatMatMultNumeric_seqaijmkl_seqdense_C",NULL);CHKERRQ(ierr); 1434 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatMatMultSymbolic_nest_seqdense_C",NULL);CHKERRQ(ierr); 1435 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatMatMultNumeric_nest_seqdense_C",NULL);CHKERRQ(ierr); 1436 1437 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatTransposeMatMultSymbolic_seqaij_seqdense_C",NULL);CHKERRQ(ierr); 1438 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatTransposeMatMultNumeric_seqaij_seqdense_C",NULL);CHKERRQ(ierr); 1439 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatTransposeMatMultSymbolic_seqaijperm_seqdense_C",NULL);CHKERRQ(ierr); 1440 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatTransposeMatMultNumeric_seqaijperm_seqdense_C",NULL);CHKERRQ(ierr); 1441 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatTransposeMatMultSymbolic_seqaijsell_seqdense_C",NULL);CHKERRQ(ierr); 1442 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatTransposeMatMultNumeric_seqaijsell_seqdense_C",NULL);CHKERRQ(ierr); 1443 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatTransposeMatMultSymbolic_seqaijmkl_seqdense_C",NULL);CHKERRQ(ierr); 1444 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatTransposeMatMultNumeric_seqaijmkl_seqdense_C",NULL);CHKERRQ(ierr); 1445 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetColumn_C",NULL);CHKERRQ(ierr); 1446 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatDenseRestoreColumn_C",NULL);CHKERRQ(ierr); 1447 PetscFunctionReturn(0); 1448 } 1449 1450 static PetscErrorCode MatTranspose_SeqDense(Mat A,MatReuse reuse,Mat *matout) 1451 { 1452 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1453 PetscErrorCode ierr; 1454 PetscInt k,j,m,n,M; 1455 PetscScalar *v,tmp; 1456 1457 PetscFunctionBegin; 1458 m = A->rmap->n; M = mat->lda; n = A->cmap->n; 1459 if (reuse == MAT_INPLACE_MATRIX && m == n) { /* in place transpose */ 1460 ierr = MatDenseGetArray(A,&v);CHKERRQ(ierr); 1461 for (j=0; j<m; j++) { 1462 for (k=0; k<j; k++) { 1463 tmp = v[j + k*M]; 1464 v[j + k*M] = v[k + j*M]; 1465 v[k + j*M] = tmp; 1466 } 1467 } 1468 ierr = MatDenseRestoreArray(A,&v);CHKERRQ(ierr); 1469 } else { /* out-of-place transpose */ 1470 Mat tmat; 1471 Mat_SeqDense *tmatd; 1472 PetscScalar *v2; 1473 PetscInt M2; 1474 1475 if (reuse != MAT_REUSE_MATRIX) { 1476 ierr = MatCreate(PetscObjectComm((PetscObject)A),&tmat);CHKERRQ(ierr); 1477 ierr = MatSetSizes(tmat,A->cmap->n,A->rmap->n,A->cmap->n,A->rmap->n);CHKERRQ(ierr); 1478 ierr = MatSetType(tmat,((PetscObject)A)->type_name);CHKERRQ(ierr); 1479 ierr = MatSeqDenseSetPreallocation(tmat,NULL);CHKERRQ(ierr); 1480 } else tmat = *matout; 1481 1482 ierr = MatDenseGetArrayRead(A,(const PetscScalar**)&v);CHKERRQ(ierr); 1483 ierr = MatDenseGetArray(tmat,&v2);CHKERRQ(ierr); 1484 tmatd = (Mat_SeqDense*)tmat->data; 1485 M2 = tmatd->lda; 1486 for (j=0; j<n; j++) { 1487 for (k=0; k<m; k++) v2[j + k*M2] = v[k + j*M]; 1488 } 1489 ierr = MatDenseRestoreArray(tmat,&v2);CHKERRQ(ierr); 1490 ierr = MatDenseRestoreArrayRead(A,(const PetscScalar**)&v);CHKERRQ(ierr); 1491 ierr = MatAssemblyBegin(tmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1492 ierr = MatAssemblyEnd(tmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1493 if (reuse == MAT_INITIAL_MATRIX || reuse == MAT_REUSE_MATRIX) *matout = tmat; 1494 else { 1495 ierr = MatHeaderMerge(A,&tmat);CHKERRQ(ierr); 1496 } 1497 } 1498 PetscFunctionReturn(0); 1499 } 1500 1501 static PetscErrorCode MatEqual_SeqDense(Mat A1,Mat A2,PetscBool *flg) 1502 { 1503 Mat_SeqDense *mat1 = (Mat_SeqDense*)A1->data; 1504 Mat_SeqDense *mat2 = (Mat_SeqDense*)A2->data; 1505 PetscInt i; 1506 const PetscScalar *v1,*v2; 1507 PetscErrorCode ierr; 1508 1509 PetscFunctionBegin; 1510 if (A1->rmap->n != A2->rmap->n) {*flg = PETSC_FALSE; PetscFunctionReturn(0);} 1511 if (A1->cmap->n != A2->cmap->n) {*flg = PETSC_FALSE; PetscFunctionReturn(0);} 1512 ierr = MatDenseGetArrayRead(A1,&v1);CHKERRQ(ierr); 1513 ierr = MatDenseGetArrayRead(A2,&v2);CHKERRQ(ierr); 1514 for (i=0; i<A1->cmap->n; i++) { 1515 ierr = PetscArraycmp(v1,v2,A1->rmap->n,flg);CHKERRQ(ierr); 1516 if (*flg == PETSC_FALSE) PetscFunctionReturn(0); 1517 v1 += mat1->lda; 1518 v2 += mat2->lda; 1519 } 1520 ierr = MatDenseRestoreArrayRead(A1,&v1);CHKERRQ(ierr); 1521 ierr = MatDenseRestoreArrayRead(A2,&v2);CHKERRQ(ierr); 1522 *flg = PETSC_TRUE; 1523 PetscFunctionReturn(0); 1524 } 1525 1526 static PetscErrorCode MatGetDiagonal_SeqDense(Mat A,Vec v) 1527 { 1528 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1529 PetscInt i,n,len; 1530 PetscScalar *x; 1531 const PetscScalar *vv; 1532 PetscErrorCode ierr; 1533 1534 PetscFunctionBegin; 1535 ierr = VecGetSize(v,&n);CHKERRQ(ierr); 1536 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 1537 len = PetscMin(A->rmap->n,A->cmap->n); 1538 ierr = MatDenseGetArrayRead(A,&vv);CHKERRQ(ierr); 1539 if (n != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming mat and vec"); 1540 for (i=0; i<len; i++) { 1541 x[i] = vv[i*mat->lda + i]; 1542 } 1543 ierr = MatDenseRestoreArrayRead(A,&vv);CHKERRQ(ierr); 1544 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 1545 PetscFunctionReturn(0); 1546 } 1547 1548 static PetscErrorCode MatDiagonalScale_SeqDense(Mat A,Vec ll,Vec rr) 1549 { 1550 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1551 const PetscScalar *l,*r; 1552 PetscScalar x,*v,*vv; 1553 PetscErrorCode ierr; 1554 PetscInt i,j,m = A->rmap->n,n = A->cmap->n; 1555 1556 PetscFunctionBegin; 1557 ierr = MatDenseGetArray(A,&vv);CHKERRQ(ierr); 1558 if (ll) { 1559 ierr = VecGetSize(ll,&m);CHKERRQ(ierr); 1560 ierr = VecGetArrayRead(ll,&l);CHKERRQ(ierr); 1561 if (m != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Left scaling vec wrong size"); 1562 for (i=0; i<m; i++) { 1563 x = l[i]; 1564 v = vv + i; 1565 for (j=0; j<n; j++) { (*v) *= x; v+= mat->lda;} 1566 } 1567 ierr = VecRestoreArrayRead(ll,&l);CHKERRQ(ierr); 1568 ierr = PetscLogFlops(1.0*n*m);CHKERRQ(ierr); 1569 } 1570 if (rr) { 1571 ierr = VecGetSize(rr,&n);CHKERRQ(ierr); 1572 ierr = VecGetArrayRead(rr,&r);CHKERRQ(ierr); 1573 if (n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Right scaling vec wrong size"); 1574 for (i=0; i<n; i++) { 1575 x = r[i]; 1576 v = vv + i*mat->lda; 1577 for (j=0; j<m; j++) (*v++) *= x; 1578 } 1579 ierr = VecRestoreArrayRead(rr,&r);CHKERRQ(ierr); 1580 ierr = PetscLogFlops(1.0*n*m);CHKERRQ(ierr); 1581 } 1582 ierr = MatDenseRestoreArray(A,&vv);CHKERRQ(ierr); 1583 PetscFunctionReturn(0); 1584 } 1585 1586 PetscErrorCode MatNorm_SeqDense(Mat A,NormType type,PetscReal *nrm) 1587 { 1588 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1589 PetscScalar *v,*vv; 1590 PetscReal sum = 0.0; 1591 PetscInt lda =mat->lda,m=A->rmap->n,i,j; 1592 PetscErrorCode ierr; 1593 1594 PetscFunctionBegin; 1595 ierr = MatDenseGetArrayRead(A,(const PetscScalar**)&vv);CHKERRQ(ierr); 1596 v = vv; 1597 if (type == NORM_FROBENIUS) { 1598 if (lda>m) { 1599 for (j=0; j<A->cmap->n; j++) { 1600 v = vv+j*lda; 1601 for (i=0; i<m; i++) { 1602 sum += PetscRealPart(PetscConj(*v)*(*v)); v++; 1603 } 1604 } 1605 } else { 1606 #if defined(PETSC_USE_REAL___FP16) 1607 PetscBLASInt one = 1,cnt = A->cmap->n*A->rmap->n; 1608 *nrm = BLASnrm2_(&cnt,v,&one); 1609 } 1610 #else 1611 for (i=0; i<A->cmap->n*A->rmap->n; i++) { 1612 sum += PetscRealPart(PetscConj(*v)*(*v)); v++; 1613 } 1614 } 1615 *nrm = PetscSqrtReal(sum); 1616 #endif 1617 ierr = PetscLogFlops(2.0*A->cmap->n*A->rmap->n);CHKERRQ(ierr); 1618 } else if (type == NORM_1) { 1619 *nrm = 0.0; 1620 for (j=0; j<A->cmap->n; j++) { 1621 v = vv + j*mat->lda; 1622 sum = 0.0; 1623 for (i=0; i<A->rmap->n; i++) { 1624 sum += PetscAbsScalar(*v); v++; 1625 } 1626 if (sum > *nrm) *nrm = sum; 1627 } 1628 ierr = PetscLogFlops(1.0*A->cmap->n*A->rmap->n);CHKERRQ(ierr); 1629 } else if (type == NORM_INFINITY) { 1630 *nrm = 0.0; 1631 for (j=0; j<A->rmap->n; j++) { 1632 v = vv + j; 1633 sum = 0.0; 1634 for (i=0; i<A->cmap->n; i++) { 1635 sum += PetscAbsScalar(*v); v += mat->lda; 1636 } 1637 if (sum > *nrm) *nrm = sum; 1638 } 1639 ierr = PetscLogFlops(1.0*A->cmap->n*A->rmap->n);CHKERRQ(ierr); 1640 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No two norm"); 1641 ierr = MatDenseRestoreArrayRead(A,(const PetscScalar**)&vv);CHKERRQ(ierr); 1642 PetscFunctionReturn(0); 1643 } 1644 1645 static PetscErrorCode MatSetOption_SeqDense(Mat A,MatOption op,PetscBool flg) 1646 { 1647 Mat_SeqDense *aij = (Mat_SeqDense*)A->data; 1648 PetscErrorCode ierr; 1649 1650 PetscFunctionBegin; 1651 switch (op) { 1652 case MAT_ROW_ORIENTED: 1653 aij->roworiented = flg; 1654 break; 1655 case MAT_NEW_NONZERO_LOCATIONS: 1656 case MAT_NEW_NONZERO_LOCATION_ERR: 1657 case MAT_NEW_NONZERO_ALLOCATION_ERR: 1658 case MAT_NEW_DIAGONALS: 1659 case MAT_KEEP_NONZERO_PATTERN: 1660 case MAT_IGNORE_OFF_PROC_ENTRIES: 1661 case MAT_USE_HASH_TABLE: 1662 case MAT_IGNORE_ZERO_ENTRIES: 1663 case MAT_IGNORE_LOWER_TRIANGULAR: 1664 case MAT_SORTED_FULL: 1665 ierr = PetscInfo1(A,"Option %s ignored\n",MatOptions[op]);CHKERRQ(ierr); 1666 break; 1667 case MAT_SPD: 1668 case MAT_SYMMETRIC: 1669 case MAT_STRUCTURALLY_SYMMETRIC: 1670 case MAT_HERMITIAN: 1671 case MAT_SYMMETRY_ETERNAL: 1672 /* These options are handled directly by MatSetOption() */ 1673 break; 1674 default: 1675 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"unknown option %s",MatOptions[op]); 1676 } 1677 PetscFunctionReturn(0); 1678 } 1679 1680 static PetscErrorCode MatZeroEntries_SeqDense(Mat A) 1681 { 1682 Mat_SeqDense *l = (Mat_SeqDense*)A->data; 1683 PetscErrorCode ierr; 1684 PetscInt lda=l->lda,m=A->rmap->n,j; 1685 PetscScalar *v; 1686 1687 PetscFunctionBegin; 1688 ierr = MatDenseGetArray(A,&v);CHKERRQ(ierr); 1689 if (lda>m) { 1690 for (j=0; j<A->cmap->n; j++) { 1691 ierr = PetscArrayzero(v+j*lda,m);CHKERRQ(ierr); 1692 } 1693 } else { 1694 ierr = PetscArrayzero(v,A->rmap->n*A->cmap->n);CHKERRQ(ierr); 1695 } 1696 ierr = MatDenseRestoreArray(A,&v);CHKERRQ(ierr); 1697 PetscFunctionReturn(0); 1698 } 1699 1700 static PetscErrorCode MatZeroRows_SeqDense(Mat A,PetscInt N,const PetscInt rows[],PetscScalar diag,Vec x,Vec b) 1701 { 1702 PetscErrorCode ierr; 1703 Mat_SeqDense *l = (Mat_SeqDense*)A->data; 1704 PetscInt m = l->lda, n = A->cmap->n, i,j; 1705 PetscScalar *slot,*bb,*v; 1706 const PetscScalar *xx; 1707 1708 PetscFunctionBegin; 1709 if (PetscDefined(USE_DEBUG)) { 1710 for (i=0; i<N; i++) { 1711 if (rows[i] < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative row requested to be zeroed"); 1712 if (rows[i] >= A->rmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row %D requested to be zeroed greater than or equal number of rows %D",rows[i],A->rmap->n); 1713 } 1714 } 1715 if (!N) PetscFunctionReturn(0); 1716 1717 /* fix right hand side if needed */ 1718 if (x && b) { 1719 ierr = VecGetArrayRead(x,&xx);CHKERRQ(ierr); 1720 ierr = VecGetArray(b,&bb);CHKERRQ(ierr); 1721 for (i=0; i<N; i++) bb[rows[i]] = diag*xx[rows[i]]; 1722 ierr = VecRestoreArrayRead(x,&xx);CHKERRQ(ierr); 1723 ierr = VecRestoreArray(b,&bb);CHKERRQ(ierr); 1724 } 1725 1726 ierr = MatDenseGetArray(A,&v);CHKERRQ(ierr); 1727 for (i=0; i<N; i++) { 1728 slot = v + rows[i]; 1729 for (j=0; j<n; j++) { *slot = 0.0; slot += m;} 1730 } 1731 if (diag != 0.0) { 1732 if (A->rmap->n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only coded for square matrices"); 1733 for (i=0; i<N; i++) { 1734 slot = v + (m+1)*rows[i]; 1735 *slot = diag; 1736 } 1737 } 1738 ierr = MatDenseRestoreArray(A,&v);CHKERRQ(ierr); 1739 PetscFunctionReturn(0); 1740 } 1741 1742 static PetscErrorCode MatDenseGetLDA_SeqDense(Mat A,PetscInt *lda) 1743 { 1744 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1745 1746 PetscFunctionBegin; 1747 *lda = mat->lda; 1748 PetscFunctionReturn(0); 1749 } 1750 1751 PetscErrorCode MatDenseGetArray_SeqDense(Mat A,PetscScalar **array) 1752 { 1753 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1754 1755 PetscFunctionBegin; 1756 *array = mat->v; 1757 PetscFunctionReturn(0); 1758 } 1759 1760 PetscErrorCode MatDenseRestoreArray_SeqDense(Mat A,PetscScalar **array) 1761 { 1762 PetscFunctionBegin; 1763 *array = NULL; 1764 PetscFunctionReturn(0); 1765 } 1766 1767 /*@C 1768 MatDenseGetLDA - gets the leading dimension of the array returned from MatDenseGetArray() 1769 1770 Logically Collective on Mat 1771 1772 Input Parameter: 1773 . mat - a MATSEQDENSE or MATMPIDENSE matrix 1774 1775 Output Parameter: 1776 . lda - the leading dimension 1777 1778 Level: intermediate 1779 1780 .seealso: MatDenseGetArray(), MatDenseRestoreArray(), MatDenseGetArrayRead(), MatDenseRestoreArrayRead(), MatSeqDenseSetLDA() 1781 @*/ 1782 PetscErrorCode MatDenseGetLDA(Mat A,PetscInt *lda) 1783 { 1784 PetscErrorCode ierr; 1785 1786 PetscFunctionBegin; 1787 ierr = PetscUseMethod(A,"MatDenseGetLDA_C",(Mat,PetscInt*),(A,lda));CHKERRQ(ierr); 1788 PetscFunctionReturn(0); 1789 } 1790 1791 /*@C 1792 MatDenseGetArray - gives access to the array where the data for a SeqDense matrix is stored 1793 1794 Logically Collective on Mat 1795 1796 Input Parameter: 1797 . mat - a MATSEQDENSE or MATMPIDENSE matrix 1798 1799 Output Parameter: 1800 . array - pointer to the data 1801 1802 Level: intermediate 1803 1804 .seealso: MatDenseRestoreArray(), MatDenseGetArrayRead(), MatDenseRestoreArrayRead() 1805 @*/ 1806 PetscErrorCode MatDenseGetArray(Mat A,PetscScalar **array) 1807 { 1808 PetscErrorCode ierr; 1809 1810 PetscFunctionBegin; 1811 ierr = PetscUseMethod(A,"MatDenseGetArray_C",(Mat,PetscScalar**),(A,array));CHKERRQ(ierr); 1812 PetscFunctionReturn(0); 1813 } 1814 1815 /*@C 1816 MatDenseRestoreArray - returns access to the array where the data for a dense matrix is stored obtained by MatDenseGetArray() 1817 1818 Logically Collective on Mat 1819 1820 Input Parameters: 1821 + mat - a MATSEQDENSE or MATMPIDENSE matrix 1822 - array - pointer to the data 1823 1824 Level: intermediate 1825 1826 .seealso: MatDenseGetArray(), MatDenseGetArrayRead(), MatDenseRestoreArrayRead() 1827 @*/ 1828 PetscErrorCode MatDenseRestoreArray(Mat A,PetscScalar **array) 1829 { 1830 PetscErrorCode ierr; 1831 1832 PetscFunctionBegin; 1833 ierr = PetscUseMethod(A,"MatDenseRestoreArray_C",(Mat,PetscScalar**),(A,array));CHKERRQ(ierr); 1834 ierr = PetscObjectStateIncrease((PetscObject)A);CHKERRQ(ierr); 1835 #if defined(PETSC_HAVE_CUDA) 1836 A->offloadmask = PETSC_OFFLOAD_CPU; 1837 #endif 1838 PetscFunctionReturn(0); 1839 } 1840 1841 /*@C 1842 MatDenseGetArrayRead - gives access to the array where the data for a SeqDense matrix is stored 1843 1844 Not Collective 1845 1846 Input Parameter: 1847 . mat - a MATSEQDENSE or MATMPIDENSE matrix 1848 1849 Output Parameter: 1850 . array - pointer to the data 1851 1852 Level: intermediate 1853 1854 .seealso: MatDenseRestoreArray(), MatDenseGetArray(), MatDenseRestoreArrayRead() 1855 @*/ 1856 PetscErrorCode MatDenseGetArrayRead(Mat A,const PetscScalar **array) 1857 { 1858 PetscErrorCode ierr; 1859 1860 PetscFunctionBegin; 1861 ierr = PetscUseMethod(A,"MatDenseGetArrayRead_C",(Mat,const PetscScalar**),(A,array));CHKERRQ(ierr); 1862 PetscFunctionReturn(0); 1863 } 1864 1865 /*@C 1866 MatDenseRestoreArrayRead - returns access to the array where the data for a dense matrix is stored obtained by MatDenseGetArray() 1867 1868 Not Collective 1869 1870 Input Parameters: 1871 + mat - a MATSEQDENSE or MATMPIDENSE matrix 1872 - array - pointer to the data 1873 1874 Level: intermediate 1875 1876 .seealso: MatDenseGetArray(), MatDenseGetArrayRead(), MatDenseRestoreArray() 1877 @*/ 1878 PetscErrorCode MatDenseRestoreArrayRead(Mat A,const PetscScalar **array) 1879 { 1880 PetscErrorCode ierr; 1881 1882 PetscFunctionBegin; 1883 ierr = PetscUseMethod(A,"MatDenseRestoreArrayRead_C",(Mat,const PetscScalar**),(A,array));CHKERRQ(ierr); 1884 PetscFunctionReturn(0); 1885 } 1886 1887 static PetscErrorCode MatCreateSubMatrix_SeqDense(Mat A,IS isrow,IS iscol,PetscInt cs,MatReuse scall,Mat *B) 1888 { 1889 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1890 PetscErrorCode ierr; 1891 PetscInt i,j,nrows,ncols,blda; 1892 const PetscInt *irow,*icol; 1893 PetscScalar *av,*bv,*v = mat->v; 1894 Mat newmat; 1895 1896 PetscFunctionBegin; 1897 ierr = ISGetIndices(isrow,&irow);CHKERRQ(ierr); 1898 ierr = ISGetIndices(iscol,&icol);CHKERRQ(ierr); 1899 ierr = ISGetLocalSize(isrow,&nrows);CHKERRQ(ierr); 1900 ierr = ISGetLocalSize(iscol,&ncols);CHKERRQ(ierr); 1901 1902 /* Check submatrixcall */ 1903 if (scall == MAT_REUSE_MATRIX) { 1904 PetscInt n_cols,n_rows; 1905 ierr = MatGetSize(*B,&n_rows,&n_cols);CHKERRQ(ierr); 1906 if (n_rows != nrows || n_cols != ncols) { 1907 /* resize the result matrix to match number of requested rows/columns */ 1908 ierr = MatSetSizes(*B,nrows,ncols,nrows,ncols);CHKERRQ(ierr); 1909 } 1910 newmat = *B; 1911 } else { 1912 /* Create and fill new matrix */ 1913 ierr = MatCreate(PetscObjectComm((PetscObject)A),&newmat);CHKERRQ(ierr); 1914 ierr = MatSetSizes(newmat,nrows,ncols,nrows,ncols);CHKERRQ(ierr); 1915 ierr = MatSetType(newmat,((PetscObject)A)->type_name);CHKERRQ(ierr); 1916 ierr = MatSeqDenseSetPreallocation(newmat,NULL);CHKERRQ(ierr); 1917 } 1918 1919 /* Now extract the data pointers and do the copy,column at a time */ 1920 ierr = MatDenseGetArray(newmat,&bv);CHKERRQ(ierr); 1921 ierr = MatDenseGetLDA(newmat,&blda);CHKERRQ(ierr); 1922 for (i=0; i<ncols; i++) { 1923 av = v + mat->lda*icol[i]; 1924 for (j=0; j<nrows; j++) bv[j] = av[irow[j]]; 1925 bv += blda; 1926 } 1927 ierr = MatDenseRestoreArray(newmat,&bv);CHKERRQ(ierr); 1928 1929 /* Assemble the matrices so that the correct flags are set */ 1930 ierr = MatAssemblyBegin(newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1931 ierr = MatAssemblyEnd(newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1932 1933 /* Free work space */ 1934 ierr = ISRestoreIndices(isrow,&irow);CHKERRQ(ierr); 1935 ierr = ISRestoreIndices(iscol,&icol);CHKERRQ(ierr); 1936 *B = newmat; 1937 PetscFunctionReturn(0); 1938 } 1939 1940 static PetscErrorCode MatCreateSubMatrices_SeqDense(Mat A,PetscInt n,const IS irow[],const IS icol[],MatReuse scall,Mat *B[]) 1941 { 1942 PetscErrorCode ierr; 1943 PetscInt i; 1944 1945 PetscFunctionBegin; 1946 if (scall == MAT_INITIAL_MATRIX) { 1947 ierr = PetscCalloc1(n+1,B);CHKERRQ(ierr); 1948 } 1949 1950 for (i=0; i<n; i++) { 1951 ierr = MatCreateSubMatrix_SeqDense(A,irow[i],icol[i],PETSC_DECIDE,scall,&(*B)[i]);CHKERRQ(ierr); 1952 } 1953 PetscFunctionReturn(0); 1954 } 1955 1956 static PetscErrorCode MatAssemblyBegin_SeqDense(Mat mat,MatAssemblyType mode) 1957 { 1958 PetscFunctionBegin; 1959 PetscFunctionReturn(0); 1960 } 1961 1962 static PetscErrorCode MatAssemblyEnd_SeqDense(Mat mat,MatAssemblyType mode) 1963 { 1964 PetscFunctionBegin; 1965 PetscFunctionReturn(0); 1966 } 1967 1968 static PetscErrorCode MatCopy_SeqDense(Mat A,Mat B,MatStructure str) 1969 { 1970 Mat_SeqDense *a = (Mat_SeqDense*)A->data,*b = (Mat_SeqDense*)B->data; 1971 PetscErrorCode ierr; 1972 const PetscScalar *va; 1973 PetscScalar *vb; 1974 PetscInt lda1=a->lda,lda2=b->lda, m=A->rmap->n,n=A->cmap->n, j; 1975 1976 PetscFunctionBegin; 1977 /* If the two matrices don't have the same copy implementation, they aren't compatible for fast copy. */ 1978 if (A->ops->copy != B->ops->copy) { 1979 ierr = MatCopy_Basic(A,B,str);CHKERRQ(ierr); 1980 PetscFunctionReturn(0); 1981 } 1982 if (m != B->rmap->n || n != B->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"size(B) != size(A)"); 1983 ierr = MatDenseGetArrayRead(A,&va);CHKERRQ(ierr); 1984 ierr = MatDenseGetArray(B,&vb);CHKERRQ(ierr); 1985 if (lda1>m || lda2>m) { 1986 for (j=0; j<n; j++) { 1987 ierr = PetscArraycpy(vb+j*lda2,va+j*lda1,m);CHKERRQ(ierr); 1988 } 1989 } else { 1990 ierr = PetscArraycpy(vb,va,A->rmap->n*A->cmap->n);CHKERRQ(ierr); 1991 } 1992 ierr = MatDenseRestoreArray(B,&vb);CHKERRQ(ierr); 1993 ierr = MatDenseRestoreArrayRead(A,&va);CHKERRQ(ierr); 1994 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1995 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1996 PetscFunctionReturn(0); 1997 } 1998 1999 static PetscErrorCode MatSetUp_SeqDense(Mat A) 2000 { 2001 PetscErrorCode ierr; 2002 2003 PetscFunctionBegin; 2004 ierr = PetscLayoutSetUp(A->rmap);CHKERRQ(ierr); 2005 ierr = PetscLayoutSetUp(A->cmap);CHKERRQ(ierr); 2006 if (!A->preallocated) { 2007 ierr = MatSeqDenseSetPreallocation(A,0);CHKERRQ(ierr); 2008 } 2009 PetscFunctionReturn(0); 2010 } 2011 2012 static PetscErrorCode MatConjugate_SeqDense(Mat A) 2013 { 2014 PetscInt i,nz = A->rmap->n*A->cmap->n; 2015 PetscScalar *aa; 2016 PetscErrorCode ierr; 2017 2018 PetscFunctionBegin; 2019 ierr = MatDenseGetArray(A,&aa);CHKERRQ(ierr); 2020 for (i=0; i<nz; i++) aa[i] = PetscConj(aa[i]); 2021 ierr = MatDenseRestoreArray(A,&aa);CHKERRQ(ierr); 2022 PetscFunctionReturn(0); 2023 } 2024 2025 static PetscErrorCode MatRealPart_SeqDense(Mat A) 2026 { 2027 PetscInt i,nz = A->rmap->n*A->cmap->n; 2028 PetscScalar *aa; 2029 PetscErrorCode ierr; 2030 2031 PetscFunctionBegin; 2032 ierr = MatDenseGetArray(A,&aa);CHKERRQ(ierr); 2033 for (i=0; i<nz; i++) aa[i] = PetscRealPart(aa[i]); 2034 ierr = MatDenseRestoreArray(A,&aa);CHKERRQ(ierr); 2035 PetscFunctionReturn(0); 2036 } 2037 2038 static PetscErrorCode MatImaginaryPart_SeqDense(Mat A) 2039 { 2040 PetscInt i,nz = A->rmap->n*A->cmap->n; 2041 PetscScalar *aa; 2042 PetscErrorCode ierr; 2043 2044 PetscFunctionBegin; 2045 ierr = MatDenseGetArray(A,&aa);CHKERRQ(ierr); 2046 for (i=0; i<nz; i++) aa[i] = PetscImaginaryPart(aa[i]); 2047 ierr = MatDenseRestoreArray(A,&aa);CHKERRQ(ierr); 2048 PetscFunctionReturn(0); 2049 } 2050 2051 /* ----------------------------------------------------------------*/ 2052 PetscErrorCode MatMatMultSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat C) 2053 { 2054 PetscErrorCode ierr; 2055 PetscInt m=A->rmap->n,n=B->cmap->n; 2056 PetscBool flg; 2057 2058 PetscFunctionBegin; 2059 ierr = MatSetSizes(C,m,n,m,n);CHKERRQ(ierr); 2060 ierr = PetscObjectTypeCompare((PetscObject)B,((PetscObject)A)->type_name,&flg);CHKERRQ(ierr); 2061 ierr = MatSetType(C,flg ? ((PetscObject)A)->type_name : MATDENSE);CHKERRQ(ierr); 2062 ierr = MatSetUp(C);CHKERRQ(ierr); 2063 PetscFunctionReturn(0); 2064 } 2065 2066 PetscErrorCode MatMatMultNumeric_SeqDense_SeqDense(Mat A,Mat B,Mat C) 2067 { 2068 Mat_SeqDense *a,*b=(Mat_SeqDense*)B->data,*c=(Mat_SeqDense*)C->data; 2069 PetscBLASInt m,n,k; 2070 const PetscScalar *av,*bv; 2071 PetscScalar *cv; 2072 PetscScalar _DOne=1.0,_DZero=0.0; 2073 PetscBool flg; 2074 PetscErrorCode (*numeric)(Mat,Mat,Mat)=NULL; 2075 PetscErrorCode ierr; 2076 2077 PetscFunctionBegin; 2078 /* Handle case where where user provided the final C matrix rather than calling MatMatMult() with MAT_INITIAL_MATRIX*/ 2079 ierr = PetscObjectTypeCompare((PetscObject)A,MATSEQAIJ,&flg);CHKERRQ(ierr); 2080 if (flg) numeric = MatMatMultNumeric_SeqAIJ_SeqDense; 2081 ierr = PetscObjectTypeCompare((PetscObject)A,MATSEQBAIJ,&flg);CHKERRQ(ierr); 2082 if (flg) numeric = MatMatMultNumeric_SeqBAIJ_SeqDense; 2083 ierr = PetscObjectTypeCompare((PetscObject)A,MATSEQSBAIJ,&flg);CHKERRQ(ierr); 2084 if (flg) numeric = MatMatMultNumeric_SeqSBAIJ_SeqDense; 2085 ierr = PetscObjectTypeCompare((PetscObject)A,MATNEST,&flg);CHKERRQ(ierr); 2086 if (flg) numeric = MatMatMultNumeric_Nest_Dense; 2087 if (numeric) { 2088 C->ops->matmultnumeric = numeric; 2089 ierr = (*numeric)(A,B,C);CHKERRQ(ierr); 2090 PetscFunctionReturn(0); 2091 } 2092 a = (Mat_SeqDense*)A->data; 2093 ierr = PetscBLASIntCast(C->rmap->n,&m);CHKERRQ(ierr); 2094 ierr = PetscBLASIntCast(C->cmap->n,&n);CHKERRQ(ierr); 2095 ierr = PetscBLASIntCast(A->cmap->n,&k);CHKERRQ(ierr); 2096 if (!m || !n || !k) PetscFunctionReturn(0); 2097 ierr = MatDenseGetArrayRead(A,&av);CHKERRQ(ierr); 2098 ierr = MatDenseGetArrayRead(B,&bv);CHKERRQ(ierr); 2099 ierr = MatDenseGetArray(C,&cv);CHKERRQ(ierr); 2100 PetscStackCallBLAS("BLASgemm",BLASgemm_("N","N",&m,&n,&k,&_DOne,av,&a->lda,bv,&b->lda,&_DZero,cv,&c->lda)); 2101 ierr = PetscLogFlops(1.0*m*n*k + 1.0*m*n*(k-1));CHKERRQ(ierr); 2102 ierr = MatDenseRestoreArrayRead(A,&av);CHKERRQ(ierr); 2103 ierr = MatDenseRestoreArrayRead(B,&bv);CHKERRQ(ierr); 2104 ierr = MatDenseRestoreArray(C,&cv);CHKERRQ(ierr); 2105 PetscFunctionReturn(0); 2106 } 2107 2108 PetscErrorCode MatMatTransposeMultSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat C) 2109 { 2110 PetscErrorCode ierr; 2111 PetscInt m=A->rmap->n,n=B->rmap->n; 2112 PetscBool flg; 2113 2114 PetscFunctionBegin; 2115 ierr = MatSetSizes(C,m,n,m,n);CHKERRQ(ierr); 2116 ierr = PetscObjectTypeCompare((PetscObject)B,((PetscObject)A)->type_name,&flg);CHKERRQ(ierr); 2117 ierr = MatSetType(C,flg ? ((PetscObject)A)->type_name : MATDENSE);CHKERRQ(ierr); 2118 ierr = MatSetUp(C);CHKERRQ(ierr); 2119 PetscFunctionReturn(0); 2120 } 2121 2122 PetscErrorCode MatMatTransposeMultNumeric_SeqDense_SeqDense(Mat A,Mat B,Mat C) 2123 { 2124 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 2125 Mat_SeqDense *b = (Mat_SeqDense*)B->data; 2126 Mat_SeqDense *c = (Mat_SeqDense*)C->data; 2127 PetscBLASInt m,n,k; 2128 PetscScalar _DOne=1.0,_DZero=0.0; 2129 PetscErrorCode ierr; 2130 2131 PetscFunctionBegin; 2132 ierr = PetscBLASIntCast(C->rmap->n,&m);CHKERRQ(ierr); 2133 ierr = PetscBLASIntCast(C->cmap->n,&n);CHKERRQ(ierr); 2134 ierr = PetscBLASIntCast(A->cmap->n,&k);CHKERRQ(ierr); 2135 if (!m || !n || !k) PetscFunctionReturn(0); 2136 PetscStackCallBLAS("BLASgemm",BLASgemm_("N","T",&m,&n,&k,&_DOne,a->v,&a->lda,b->v,&b->lda,&_DZero,c->v,&c->lda)); 2137 ierr = PetscLogFlops(1.0*m*n*k + 1.0*m*n*(k-1));CHKERRQ(ierr); 2138 PetscFunctionReturn(0); 2139 } 2140 2141 PetscErrorCode MatTransposeMatMultSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat C) 2142 { 2143 PetscErrorCode ierr; 2144 PetscInt m=A->cmap->n,n=B->cmap->n; 2145 PetscBool flg; 2146 2147 PetscFunctionBegin; 2148 ierr = MatSetSizes(C,m,n,m,n);CHKERRQ(ierr); 2149 ierr = PetscObjectTypeCompare((PetscObject)B,((PetscObject)A)->type_name,&flg);CHKERRQ(ierr); 2150 ierr = MatSetType(C,flg ? ((PetscObject)A)->type_name : MATDENSE);CHKERRQ(ierr); 2151 ierr = MatSetUp(C);CHKERRQ(ierr); 2152 PetscFunctionReturn(0); 2153 } 2154 2155 PetscErrorCode MatTransposeMatMultNumeric_SeqDense_SeqDense(Mat A,Mat B,Mat C) 2156 { 2157 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 2158 Mat_SeqDense *b = (Mat_SeqDense*)B->data; 2159 Mat_SeqDense *c = (Mat_SeqDense*)C->data; 2160 PetscBLASInt m,n,k; 2161 PetscScalar _DOne=1.0,_DZero=0.0; 2162 PetscErrorCode ierr; 2163 2164 PetscFunctionBegin; 2165 ierr = PetscBLASIntCast(C->rmap->n,&m);CHKERRQ(ierr); 2166 ierr = PetscBLASIntCast(C->cmap->n,&n);CHKERRQ(ierr); 2167 ierr = PetscBLASIntCast(A->rmap->n,&k);CHKERRQ(ierr); 2168 if (!m || !n || !k) PetscFunctionReturn(0); 2169 PetscStackCallBLAS("BLASgemm",BLASgemm_("T","N",&m,&n,&k,&_DOne,a->v,&a->lda,b->v,&b->lda,&_DZero,c->v,&c->lda)); 2170 ierr = PetscLogFlops(1.0*m*n*k + 1.0*m*n*(k-1));CHKERRQ(ierr); 2171 PetscFunctionReturn(0); 2172 } 2173 2174 /* ----------------------------------------------- */ 2175 static PetscErrorCode MatProductSetFromOptions_SeqDense_AB(Mat C) 2176 { 2177 PetscFunctionBegin; 2178 C->ops->matmultsymbolic = MatMatMultSymbolic_SeqDense_SeqDense; 2179 C->ops->productsymbolic = MatProductSymbolic_AB; 2180 /* dense mat may not call MatProductSymbolic(), thus set C->ops->productnumeric here */ 2181 C->ops->productnumeric = MatProductNumeric_AB; 2182 PetscFunctionReturn(0); 2183 } 2184 2185 static PetscErrorCode MatProductSetFromOptions_SeqDense_AtB(Mat C) 2186 { 2187 PetscFunctionBegin; 2188 C->ops->transposematmultsymbolic = MatTransposeMatMultSymbolic_SeqDense_SeqDense; 2189 C->ops->productsymbolic = MatProductSymbolic_AtB; 2190 C->ops->productnumeric = MatProductNumeric_AtB; 2191 PetscFunctionReturn(0); 2192 } 2193 2194 static PetscErrorCode MatProductSetFromOptions_SeqDense_ABt(Mat C) 2195 { 2196 PetscFunctionBegin; 2197 C->ops->mattransposemultsymbolic = MatMatTransposeMultSymbolic_SeqDense_SeqDense; 2198 C->ops->productsymbolic = MatProductSymbolic_ABt; 2199 C->ops->productnumeric = MatProductNumeric_ABt; 2200 PetscFunctionReturn(0); 2201 } 2202 2203 static PetscErrorCode MatProductSetFromOptions_SeqDense_PtAP(Mat C) 2204 { 2205 PetscFunctionBegin; 2206 C->ops->productsymbolic = MatProductSymbolic_Basic; 2207 PetscFunctionReturn(0); 2208 } 2209 2210 PETSC_INTERN PetscErrorCode MatProductSetFromOptions_SeqDense(Mat C) 2211 { 2212 PetscErrorCode ierr; 2213 Mat_Product *product = C->product; 2214 2215 PetscFunctionBegin; 2216 switch (product->type) { 2217 case MATPRODUCT_AB: 2218 ierr = MatProductSetFromOptions_SeqDense_AB(C);CHKERRQ(ierr); 2219 break; 2220 case MATPRODUCT_AtB: 2221 ierr = MatProductSetFromOptions_SeqDense_AtB(C);CHKERRQ(ierr); 2222 break; 2223 case MATPRODUCT_ABt: 2224 ierr = MatProductSetFromOptions_SeqDense_ABt(C);CHKERRQ(ierr); 2225 break; 2226 case MATPRODUCT_PtAP: 2227 ierr = MatProductSetFromOptions_SeqDense_PtAP(C);CHKERRQ(ierr); 2228 break; 2229 default: SETERRQ1(PetscObjectComm((PetscObject)C),PETSC_ERR_SUP,"MatProduct type %s is not supported for SeqDense and SeqDense matrices",MatProductTypes[product->type]); 2230 } 2231 PetscFunctionReturn(0); 2232 } 2233 /* ----------------------------------------------- */ 2234 2235 static PetscErrorCode MatGetRowMax_SeqDense(Mat A,Vec v,PetscInt idx[]) 2236 { 2237 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 2238 PetscErrorCode ierr; 2239 PetscInt i,j,m = A->rmap->n,n = A->cmap->n,p; 2240 PetscScalar *x; 2241 const PetscScalar *aa; 2242 2243 PetscFunctionBegin; 2244 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2245 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2246 ierr = VecGetLocalSize(v,&p);CHKERRQ(ierr); 2247 ierr = MatDenseGetArrayRead(A,&aa);CHKERRQ(ierr); 2248 if (p != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 2249 for (i=0; i<m; i++) { 2250 x[i] = aa[i]; if (idx) idx[i] = 0; 2251 for (j=1; j<n; j++) { 2252 if (PetscRealPart(x[i]) < PetscRealPart(aa[i+a->lda*j])) {x[i] = aa[i + a->lda*j]; if (idx) idx[i] = j;} 2253 } 2254 } 2255 ierr = MatDenseRestoreArrayRead(A,&aa);CHKERRQ(ierr); 2256 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2257 PetscFunctionReturn(0); 2258 } 2259 2260 static PetscErrorCode MatGetRowMaxAbs_SeqDense(Mat A,Vec v,PetscInt idx[]) 2261 { 2262 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 2263 PetscErrorCode ierr; 2264 PetscInt i,j,m = A->rmap->n,n = A->cmap->n,p; 2265 PetscScalar *x; 2266 PetscReal atmp; 2267 const PetscScalar *aa; 2268 2269 PetscFunctionBegin; 2270 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2271 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2272 ierr = VecGetLocalSize(v,&p);CHKERRQ(ierr); 2273 ierr = MatDenseGetArrayRead(A,&aa);CHKERRQ(ierr); 2274 if (p != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 2275 for (i=0; i<m; i++) { 2276 x[i] = PetscAbsScalar(aa[i]); 2277 for (j=1; j<n; j++) { 2278 atmp = PetscAbsScalar(aa[i+a->lda*j]); 2279 if (PetscAbsScalar(x[i]) < atmp) {x[i] = atmp; if (idx) idx[i] = j;} 2280 } 2281 } 2282 ierr = MatDenseRestoreArrayRead(A,&aa);CHKERRQ(ierr); 2283 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2284 PetscFunctionReturn(0); 2285 } 2286 2287 static PetscErrorCode MatGetRowMin_SeqDense(Mat A,Vec v,PetscInt idx[]) 2288 { 2289 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 2290 PetscErrorCode ierr; 2291 PetscInt i,j,m = A->rmap->n,n = A->cmap->n,p; 2292 PetscScalar *x; 2293 const PetscScalar *aa; 2294 2295 PetscFunctionBegin; 2296 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2297 ierr = MatDenseGetArrayRead(A,&aa);CHKERRQ(ierr); 2298 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2299 ierr = VecGetLocalSize(v,&p);CHKERRQ(ierr); 2300 if (p != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 2301 for (i=0; i<m; i++) { 2302 x[i] = aa[i]; if (idx) idx[i] = 0; 2303 for (j=1; j<n; j++) { 2304 if (PetscRealPart(x[i]) > PetscRealPart(aa[i+a->lda*j])) {x[i] = aa[i + a->lda*j]; if (idx) idx[i] = j;} 2305 } 2306 } 2307 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2308 ierr = MatDenseRestoreArrayRead(A,&aa);CHKERRQ(ierr); 2309 PetscFunctionReturn(0); 2310 } 2311 2312 PetscErrorCode MatGetColumnVector_SeqDense(Mat A,Vec v,PetscInt col) 2313 { 2314 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 2315 PetscErrorCode ierr; 2316 PetscScalar *x; 2317 const PetscScalar *aa; 2318 2319 PetscFunctionBegin; 2320 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2321 ierr = MatDenseGetArrayRead(A,&aa);CHKERRQ(ierr); 2322 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2323 ierr = PetscArraycpy(x,aa+col*a->lda,A->rmap->n);CHKERRQ(ierr); 2324 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2325 ierr = MatDenseRestoreArrayRead(A,&aa);CHKERRQ(ierr); 2326 PetscFunctionReturn(0); 2327 } 2328 2329 PETSC_INTERN PetscErrorCode MatGetColumnNorms_SeqDense(Mat A,NormType type,PetscReal *norms) 2330 { 2331 PetscErrorCode ierr; 2332 PetscInt i,j,m,n; 2333 const PetscScalar *a; 2334 2335 PetscFunctionBegin; 2336 ierr = MatGetSize(A,&m,&n);CHKERRQ(ierr); 2337 ierr = PetscArrayzero(norms,n);CHKERRQ(ierr); 2338 ierr = MatDenseGetArrayRead(A,&a);CHKERRQ(ierr); 2339 if (type == NORM_2) { 2340 for (i=0; i<n; i++) { 2341 for (j=0; j<m; j++) { 2342 norms[i] += PetscAbsScalar(a[j]*a[j]); 2343 } 2344 a += m; 2345 } 2346 } else if (type == NORM_1) { 2347 for (i=0; i<n; i++) { 2348 for (j=0; j<m; j++) { 2349 norms[i] += PetscAbsScalar(a[j]); 2350 } 2351 a += m; 2352 } 2353 } else if (type == NORM_INFINITY) { 2354 for (i=0; i<n; i++) { 2355 for (j=0; j<m; j++) { 2356 norms[i] = PetscMax(PetscAbsScalar(a[j]),norms[i]); 2357 } 2358 a += m; 2359 } 2360 } else SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Unknown NormType"); 2361 ierr = MatDenseRestoreArrayRead(A,&a);CHKERRQ(ierr); 2362 if (type == NORM_2) { 2363 for (i=0; i<n; i++) norms[i] = PetscSqrtReal(norms[i]); 2364 } 2365 PetscFunctionReturn(0); 2366 } 2367 2368 static PetscErrorCode MatSetRandom_SeqDense(Mat x,PetscRandom rctx) 2369 { 2370 PetscErrorCode ierr; 2371 PetscScalar *a; 2372 PetscInt lda,m,n,i,j; 2373 2374 PetscFunctionBegin; 2375 ierr = MatGetSize(x,&m,&n);CHKERRQ(ierr); 2376 ierr = MatDenseGetLDA(x,&lda);CHKERRQ(ierr); 2377 ierr = MatDenseGetArray(x,&a);CHKERRQ(ierr); 2378 for (j=0; j<n; j++) { 2379 for (i=0; i<m; i++) { 2380 ierr = PetscRandomGetValue(rctx,a+j*lda+i);CHKERRQ(ierr); 2381 } 2382 } 2383 ierr = MatDenseRestoreArray(x,&a);CHKERRQ(ierr); 2384 PetscFunctionReturn(0); 2385 } 2386 2387 static PetscErrorCode MatMissingDiagonal_SeqDense(Mat A,PetscBool *missing,PetscInt *d) 2388 { 2389 PetscFunctionBegin; 2390 *missing = PETSC_FALSE; 2391 PetscFunctionReturn(0); 2392 } 2393 2394 /* vals is not const */ 2395 static PetscErrorCode MatDenseGetColumn_SeqDense(Mat A,PetscInt col,PetscScalar **vals) 2396 { 2397 PetscErrorCode ierr; 2398 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 2399 PetscScalar *v; 2400 2401 PetscFunctionBegin; 2402 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2403 ierr = MatDenseGetArray(A,&v);CHKERRQ(ierr); 2404 *vals = v+col*a->lda; 2405 ierr = MatDenseRestoreArray(A,&v);CHKERRQ(ierr); 2406 PetscFunctionReturn(0); 2407 } 2408 2409 static PetscErrorCode MatDenseRestoreColumn_SeqDense(Mat A,PetscScalar **vals) 2410 { 2411 PetscFunctionBegin; 2412 *vals = 0; /* user cannot accidently use the array later */ 2413 PetscFunctionReturn(0); 2414 } 2415 2416 /* -------------------------------------------------------------------*/ 2417 static struct _MatOps MatOps_Values = { MatSetValues_SeqDense, 2418 MatGetRow_SeqDense, 2419 MatRestoreRow_SeqDense, 2420 MatMult_SeqDense, 2421 /* 4*/ MatMultAdd_SeqDense, 2422 MatMultTranspose_SeqDense, 2423 MatMultTransposeAdd_SeqDense, 2424 0, 2425 0, 2426 0, 2427 /* 10*/ 0, 2428 MatLUFactor_SeqDense, 2429 MatCholeskyFactor_SeqDense, 2430 MatSOR_SeqDense, 2431 MatTranspose_SeqDense, 2432 /* 15*/ MatGetInfo_SeqDense, 2433 MatEqual_SeqDense, 2434 MatGetDiagonal_SeqDense, 2435 MatDiagonalScale_SeqDense, 2436 MatNorm_SeqDense, 2437 /* 20*/ MatAssemblyBegin_SeqDense, 2438 MatAssemblyEnd_SeqDense, 2439 MatSetOption_SeqDense, 2440 MatZeroEntries_SeqDense, 2441 /* 24*/ MatZeroRows_SeqDense, 2442 0, 2443 0, 2444 0, 2445 0, 2446 /* 29*/ MatSetUp_SeqDense, 2447 0, 2448 0, 2449 0, 2450 0, 2451 /* 34*/ MatDuplicate_SeqDense, 2452 0, 2453 0, 2454 0, 2455 0, 2456 /* 39*/ MatAXPY_SeqDense, 2457 MatCreateSubMatrices_SeqDense, 2458 0, 2459 MatGetValues_SeqDense, 2460 MatCopy_SeqDense, 2461 /* 44*/ MatGetRowMax_SeqDense, 2462 MatScale_SeqDense, 2463 MatShift_Basic, 2464 0, 2465 MatZeroRowsColumns_SeqDense, 2466 /* 49*/ MatSetRandom_SeqDense, 2467 0, 2468 0, 2469 0, 2470 0, 2471 /* 54*/ 0, 2472 0, 2473 0, 2474 0, 2475 0, 2476 /* 59*/ 0, 2477 MatDestroy_SeqDense, 2478 MatView_SeqDense, 2479 0, 2480 0, 2481 /* 64*/ 0, 2482 0, 2483 0, 2484 0, 2485 0, 2486 /* 69*/ MatGetRowMaxAbs_SeqDense, 2487 0, 2488 0, 2489 0, 2490 0, 2491 /* 74*/ 0, 2492 0, 2493 0, 2494 0, 2495 0, 2496 /* 79*/ 0, 2497 0, 2498 0, 2499 0, 2500 /* 83*/ MatLoad_SeqDense, 2501 MatIsSymmetric_SeqDense, 2502 MatIsHermitian_SeqDense, 2503 0, 2504 0, 2505 0, 2506 /* 89*/ 0, 2507 0, 2508 MatMatMultNumeric_SeqDense_SeqDense, 2509 0, 2510 0, 2511 /* 94*/ 0, 2512 0, 2513 0, 2514 MatMatTransposeMultNumeric_SeqDense_SeqDense, 2515 0, 2516 /* 99*/ MatProductSetFromOptions_SeqDense, 2517 0, 2518 0, 2519 MatConjugate_SeqDense, 2520 0, 2521 /*104*/ 0, 2522 MatRealPart_SeqDense, 2523 MatImaginaryPart_SeqDense, 2524 0, 2525 0, 2526 /*109*/ 0, 2527 0, 2528 MatGetRowMin_SeqDense, 2529 MatGetColumnVector_SeqDense, 2530 MatMissingDiagonal_SeqDense, 2531 /*114*/ 0, 2532 0, 2533 0, 2534 0, 2535 0, 2536 /*119*/ 0, 2537 0, 2538 0, 2539 0, 2540 0, 2541 /*124*/ 0, 2542 MatGetColumnNorms_SeqDense, 2543 0, 2544 0, 2545 0, 2546 /*129*/ 0, 2547 0, 2548 0, 2549 MatTransposeMatMultNumeric_SeqDense_SeqDense, 2550 0, 2551 /*134*/ 0, 2552 0, 2553 0, 2554 0, 2555 0, 2556 /*139*/ 0, 2557 0, 2558 0, 2559 0, 2560 0, 2561 MatCreateMPIMatConcatenateSeqMat_SeqDense, 2562 /*145*/ 0, 2563 0, 2564 0 2565 }; 2566 2567 /*@C 2568 MatCreateSeqDense - Creates a sequential dense matrix that 2569 is stored in column major order (the usual Fortran 77 manner). Many 2570 of the matrix operations use the BLAS and LAPACK routines. 2571 2572 Collective 2573 2574 Input Parameters: 2575 + comm - MPI communicator, set to PETSC_COMM_SELF 2576 . m - number of rows 2577 . n - number of columns 2578 - data - optional location of matrix data in column major order. Set data=NULL for PETSc 2579 to control all matrix memory allocation. 2580 2581 Output Parameter: 2582 . A - the matrix 2583 2584 Notes: 2585 The data input variable is intended primarily for Fortran programmers 2586 who wish to allocate their own matrix memory space. Most users should 2587 set data=NULL. 2588 2589 Level: intermediate 2590 2591 .seealso: MatCreate(), MatCreateDense(), MatSetValues() 2592 @*/ 2593 PetscErrorCode MatCreateSeqDense(MPI_Comm comm,PetscInt m,PetscInt n,PetscScalar *data,Mat *A) 2594 { 2595 PetscErrorCode ierr; 2596 2597 PetscFunctionBegin; 2598 ierr = MatCreate(comm,A);CHKERRQ(ierr); 2599 ierr = MatSetSizes(*A,m,n,m,n);CHKERRQ(ierr); 2600 ierr = MatSetType(*A,MATSEQDENSE);CHKERRQ(ierr); 2601 ierr = MatSeqDenseSetPreallocation(*A,data);CHKERRQ(ierr); 2602 PetscFunctionReturn(0); 2603 } 2604 2605 /*@C 2606 MatSeqDenseSetPreallocation - Sets the array used for storing the matrix elements 2607 2608 Collective 2609 2610 Input Parameters: 2611 + B - the matrix 2612 - data - the array (or NULL) 2613 2614 Notes: 2615 The data input variable is intended primarily for Fortran programmers 2616 who wish to allocate their own matrix memory space. Most users should 2617 need not call this routine. 2618 2619 Level: intermediate 2620 2621 .seealso: MatCreate(), MatCreateDense(), MatSetValues(), MatSeqDenseSetLDA() 2622 2623 @*/ 2624 PetscErrorCode MatSeqDenseSetPreallocation(Mat B,PetscScalar data[]) 2625 { 2626 PetscErrorCode ierr; 2627 2628 PetscFunctionBegin; 2629 ierr = PetscTryMethod(B,"MatSeqDenseSetPreallocation_C",(Mat,PetscScalar[]),(B,data));CHKERRQ(ierr); 2630 PetscFunctionReturn(0); 2631 } 2632 2633 PetscErrorCode MatSeqDenseSetPreallocation_SeqDense(Mat B,PetscScalar *data) 2634 { 2635 Mat_SeqDense *b; 2636 PetscErrorCode ierr; 2637 2638 PetscFunctionBegin; 2639 B->preallocated = PETSC_TRUE; 2640 2641 ierr = PetscLayoutSetUp(B->rmap);CHKERRQ(ierr); 2642 ierr = PetscLayoutSetUp(B->cmap);CHKERRQ(ierr); 2643 2644 b = (Mat_SeqDense*)B->data; 2645 b->Mmax = B->rmap->n; 2646 b->Nmax = B->cmap->n; 2647 if (b->lda <= 0 || b->changelda) b->lda = B->rmap->n; 2648 2649 ierr = PetscIntMultError(b->lda,b->Nmax,NULL);CHKERRQ(ierr); 2650 if (!data) { /* petsc-allocated storage */ 2651 if (!b->user_alloc) { ierr = PetscFree(b->v);CHKERRQ(ierr); } 2652 ierr = PetscCalloc1((size_t)b->lda*b->Nmax,&b->v);CHKERRQ(ierr); 2653 ierr = PetscLogObjectMemory((PetscObject)B,b->lda*b->Nmax*sizeof(PetscScalar));CHKERRQ(ierr); 2654 2655 b->user_alloc = PETSC_FALSE; 2656 } else { /* user-allocated storage */ 2657 if (!b->user_alloc) { ierr = PetscFree(b->v);CHKERRQ(ierr); } 2658 b->v = data; 2659 b->user_alloc = PETSC_TRUE; 2660 } 2661 B->assembled = PETSC_TRUE; 2662 PetscFunctionReturn(0); 2663 } 2664 2665 #if defined(PETSC_HAVE_ELEMENTAL) 2666 PETSC_INTERN PetscErrorCode MatConvert_SeqDense_Elemental(Mat A, MatType newtype,MatReuse reuse,Mat *newmat) 2667 { 2668 Mat mat_elemental; 2669 PetscErrorCode ierr; 2670 const PetscScalar *array; 2671 PetscScalar *v_colwise; 2672 PetscInt M=A->rmap->N,N=A->cmap->N,i,j,k,*rows,*cols; 2673 2674 PetscFunctionBegin; 2675 ierr = PetscMalloc3(M*N,&v_colwise,M,&rows,N,&cols);CHKERRQ(ierr); 2676 ierr = MatDenseGetArrayRead(A,&array);CHKERRQ(ierr); 2677 /* convert column-wise array into row-wise v_colwise, see MatSetValues_Elemental() */ 2678 k = 0; 2679 for (j=0; j<N; j++) { 2680 cols[j] = j; 2681 for (i=0; i<M; i++) { 2682 v_colwise[j*M+i] = array[k++]; 2683 } 2684 } 2685 for (i=0; i<M; i++) { 2686 rows[i] = i; 2687 } 2688 ierr = MatDenseRestoreArrayRead(A,&array);CHKERRQ(ierr); 2689 2690 ierr = MatCreate(PetscObjectComm((PetscObject)A), &mat_elemental);CHKERRQ(ierr); 2691 ierr = MatSetSizes(mat_elemental,PETSC_DECIDE,PETSC_DECIDE,M,N);CHKERRQ(ierr); 2692 ierr = MatSetType(mat_elemental,MATELEMENTAL);CHKERRQ(ierr); 2693 ierr = MatSetUp(mat_elemental);CHKERRQ(ierr); 2694 2695 /* PETSc-Elemental interaface uses axpy for setting off-processor entries, only ADD_VALUES is allowed */ 2696 ierr = MatSetValues(mat_elemental,M,rows,N,cols,v_colwise,ADD_VALUES);CHKERRQ(ierr); 2697 ierr = MatAssemblyBegin(mat_elemental, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2698 ierr = MatAssemblyEnd(mat_elemental, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2699 ierr = PetscFree3(v_colwise,rows,cols);CHKERRQ(ierr); 2700 2701 if (reuse == MAT_INPLACE_MATRIX) { 2702 ierr = MatHeaderReplace(A,&mat_elemental);CHKERRQ(ierr); 2703 } else { 2704 *newmat = mat_elemental; 2705 } 2706 PetscFunctionReturn(0); 2707 } 2708 #endif 2709 2710 /*@C 2711 MatSeqDenseSetLDA - Declare the leading dimension of the user-provided array 2712 2713 Input parameter: 2714 + A - the matrix 2715 - lda - the leading dimension 2716 2717 Notes: 2718 This routine is to be used in conjunction with MatSeqDenseSetPreallocation(); 2719 it asserts that the preallocation has a leading dimension (the LDA parameter 2720 of Blas and Lapack fame) larger than M, the first dimension of the matrix. 2721 2722 Level: intermediate 2723 2724 .seealso: MatCreate(), MatCreateSeqDense(), MatSeqDenseSetPreallocation(), MatSetMaximumSize() 2725 2726 @*/ 2727 PetscErrorCode MatSeqDenseSetLDA(Mat B,PetscInt lda) 2728 { 2729 Mat_SeqDense *b = (Mat_SeqDense*)B->data; 2730 2731 PetscFunctionBegin; 2732 if (lda < B->rmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"LDA %D must be at least matrix dimension %D",lda,B->rmap->n); 2733 b->lda = lda; 2734 b->changelda = PETSC_FALSE; 2735 b->Mmax = PetscMax(b->Mmax,lda); 2736 PetscFunctionReturn(0); 2737 } 2738 2739 PetscErrorCode MatCreateMPIMatConcatenateSeqMat_SeqDense(MPI_Comm comm,Mat inmat,PetscInt n,MatReuse scall,Mat *outmat) 2740 { 2741 PetscErrorCode ierr; 2742 PetscMPIInt size; 2743 2744 PetscFunctionBegin; 2745 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 2746 if (size == 1) { 2747 if (scall == MAT_INITIAL_MATRIX) { 2748 ierr = MatDuplicate(inmat,MAT_COPY_VALUES,outmat);CHKERRQ(ierr); 2749 } else { 2750 ierr = MatCopy(inmat,*outmat,SAME_NONZERO_PATTERN);CHKERRQ(ierr); 2751 } 2752 } else { 2753 ierr = MatCreateMPIMatConcatenateSeqMat_MPIDense(comm,inmat,n,scall,outmat);CHKERRQ(ierr); 2754 } 2755 PetscFunctionReturn(0); 2756 } 2757 2758 /*MC 2759 MATSEQDENSE - MATSEQDENSE = "seqdense" - A matrix type to be used for sequential dense matrices. 2760 2761 Options Database Keys: 2762 . -mat_type seqdense - sets the matrix type to "seqdense" during a call to MatSetFromOptions() 2763 2764 Level: beginner 2765 2766 .seealso: MatCreateSeqDense() 2767 2768 M*/ 2769 PetscErrorCode MatCreate_SeqDense(Mat B) 2770 { 2771 Mat_SeqDense *b; 2772 PetscErrorCode ierr; 2773 PetscMPIInt size; 2774 2775 PetscFunctionBegin; 2776 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)B),&size);CHKERRQ(ierr); 2777 if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Comm must be of size 1"); 2778 2779 ierr = PetscNewLog(B,&b);CHKERRQ(ierr); 2780 ierr = PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 2781 B->data = (void*)b; 2782 2783 b->roworiented = PETSC_TRUE; 2784 2785 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseGetLDA_C",MatDenseGetLDA_SeqDense);CHKERRQ(ierr); 2786 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseGetArray_C",MatDenseGetArray_SeqDense);CHKERRQ(ierr); 2787 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreArray_C",MatDenseRestoreArray_SeqDense);CHKERRQ(ierr); 2788 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDensePlaceArray_C",MatDensePlaceArray_SeqDense);CHKERRQ(ierr); 2789 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseResetArray_C",MatDenseResetArray_SeqDense);CHKERRQ(ierr); 2790 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseGetArrayRead_C",MatDenseGetArray_SeqDense);CHKERRQ(ierr); 2791 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreArrayRead_C",MatDenseRestoreArray_SeqDense);CHKERRQ(ierr); 2792 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqdense_seqaij_C",MatConvert_SeqDense_SeqAIJ);CHKERRQ(ierr); 2793 #if defined(PETSC_HAVE_ELEMENTAL) 2794 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqdense_elemental_C",MatConvert_SeqDense_Elemental);CHKERRQ(ierr); 2795 #endif 2796 #if defined(PETSC_HAVE_CUDA) 2797 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqdense_seqdensecuda_C",MatConvert_SeqDense_SeqDenseCUDA);CHKERRQ(ierr); 2798 ierr = PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqdensecuda_seqdensecuda_C",MatProductSetFromOptions_SeqDense);CHKERRQ(ierr); 2799 ierr = PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqdensecuda_seqdense_C",MatProductSetFromOptions_SeqDense);CHKERRQ(ierr); 2800 ierr = PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqdense_seqdensecuda_C",MatProductSetFromOptions_SeqDense);CHKERRQ(ierr); 2801 ierr = PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqaijcusparse_seqdense_C",MatProductSetFromOptions_SeqAIJ_SeqDense);CHKERRQ(ierr); 2802 #endif 2803 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqDenseSetPreallocation_C",MatSeqDenseSetPreallocation_SeqDense);CHKERRQ(ierr); 2804 ierr = PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqaij_seqdense_C",MatProductSetFromOptions_SeqAIJ_SeqDense);CHKERRQ(ierr); 2805 ierr = PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqdense_seqdense_C",MatProductSetFromOptions_SeqDense);CHKERRQ(ierr); 2806 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMultSymbolic_seqaij_seqdense_C",MatMatMultSymbolic_SeqAIJ_SeqDense);CHKERRQ(ierr); 2807 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMultNumeric_seqaij_seqdense_C",MatMatMultNumeric_SeqAIJ_SeqDense);CHKERRQ(ierr); 2808 ierr = PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqbaij_seqdense_C",MatProductSetFromOptions_SeqXBAIJ_SeqDense);CHKERRQ(ierr); 2809 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMultSymbolic_seqbaij_seqdense_C",MatMatMultSymbolic_SeqBAIJ_SeqDense);CHKERRQ(ierr); 2810 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMultNumeric_seqbaij_seqdense_C",MatMatMultNumeric_SeqBAIJ_SeqDense);CHKERRQ(ierr); 2811 ierr = PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqsbaij_seqdense_C",MatProductSetFromOptions_SeqXBAIJ_SeqDense);CHKERRQ(ierr); 2812 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMultSymbolic_seqsbaij_seqdense_C",MatMatMultSymbolic_SeqSBAIJ_SeqDense);CHKERRQ(ierr); 2813 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMultNumeric_seqsbaij_seqdense_C",MatMatMultNumeric_SeqSBAIJ_SeqDense);CHKERRQ(ierr); 2814 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMultSymbolic_seqaijperm_seqdense_C",MatMatMultSymbolic_SeqAIJ_SeqDense);CHKERRQ(ierr); 2815 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMultNumeric_seqaijperm_seqdense_C",MatMatMultNumeric_SeqAIJ_SeqDense);CHKERRQ(ierr); 2816 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMultSymbolic_seqaijsell_seqdense_C",MatMatMultSymbolic_SeqAIJ_SeqDense);CHKERRQ(ierr); 2817 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMultNumeric_seqaijsell_seqdense_C",MatMatMultNumeric_SeqAIJ_SeqDense);CHKERRQ(ierr); 2818 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMultSymbolic_seqaijmkl_seqdense_C",MatMatMultSymbolic_SeqAIJ_SeqDense);CHKERRQ(ierr); 2819 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMultNumeric_seqaijmkl_seqdense_C",MatMatMultNumeric_SeqAIJ_SeqDense);CHKERRQ(ierr); 2820 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMultSymbolic_nest_seqdense_C",MatMatMultSymbolic_Nest_Dense);CHKERRQ(ierr); 2821 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMultNumeric_nest_seqdense_C",MatMatMultNumeric_Nest_Dense);CHKERRQ(ierr); 2822 2823 ierr = PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultSymbolic_seqaij_seqdense_C",MatTransposeMatMultSymbolic_SeqAIJ_SeqDense);CHKERRQ(ierr); 2824 ierr = PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultNumeric_seqaij_seqdense_C",MatTransposeMatMultNumeric_SeqAIJ_SeqDense);CHKERRQ(ierr); 2825 ierr = PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultSymbolic_seqaijperm_seqdense_C",MatTransposeMatMultSymbolic_SeqAIJ_SeqDense);CHKERRQ(ierr); 2826 ierr = PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultNumeric_seqaijperm_seqdense_C",MatTransposeMatMultNumeric_SeqAIJ_SeqDense);CHKERRQ(ierr); 2827 ierr = PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultSymbolic_seqaijsell_seqdense_C",MatTransposeMatMultSymbolic_SeqAIJ_SeqDense);CHKERRQ(ierr); 2828 ierr = PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultNumeric_seqaijsell_seqdense_C",MatTransposeMatMultNumeric_SeqAIJ_SeqDense);CHKERRQ(ierr); 2829 2830 ierr = PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultSymbolic_seqaijmkl_seqdense_C",MatTransposeMatMultSymbolic_SeqAIJ_SeqDense);CHKERRQ(ierr); 2831 ierr = PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultNumeric_seqaijmkl_seqdense_C",MatTransposeMatMultNumeric_SeqAIJ_SeqDense);CHKERRQ(ierr); 2832 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseGetColumn_C",MatDenseGetColumn_SeqDense);CHKERRQ(ierr); 2833 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreColumn_C",MatDenseRestoreColumn_SeqDense);CHKERRQ(ierr); 2834 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQDENSE);CHKERRQ(ierr); 2835 PetscFunctionReturn(0); 2836 } 2837 2838 /*@C 2839 MatDenseGetColumn - gives access to a column of a dense matrix. This is only the local part of the column. You MUST call MatDenseRestoreColumn() to avoid memory bleeding. 2840 2841 Not Collective 2842 2843 Input Parameter: 2844 + mat - a MATSEQDENSE or MATMPIDENSE matrix 2845 - col - column index 2846 2847 Output Parameter: 2848 . vals - pointer to the data 2849 2850 Level: intermediate 2851 2852 .seealso: MatDenseRestoreColumn() 2853 @*/ 2854 PetscErrorCode MatDenseGetColumn(Mat A,PetscInt col,PetscScalar **vals) 2855 { 2856 PetscErrorCode ierr; 2857 2858 PetscFunctionBegin; 2859 ierr = PetscUseMethod(A,"MatDenseGetColumn_C",(Mat,PetscInt,PetscScalar**),(A,col,vals));CHKERRQ(ierr); 2860 PetscFunctionReturn(0); 2861 } 2862 2863 /*@C 2864 MatDenseRestoreColumn - returns access to a column of a dense matrix which is returned by MatDenseGetColumn(). 2865 2866 Not Collective 2867 2868 Input Parameter: 2869 . mat - a MATSEQDENSE or MATMPIDENSE matrix 2870 2871 Output Parameter: 2872 . vals - pointer to the data 2873 2874 Level: intermediate 2875 2876 .seealso: MatDenseGetColumn() 2877 @*/ 2878 PetscErrorCode MatDenseRestoreColumn(Mat A,PetscScalar **vals) 2879 { 2880 PetscErrorCode ierr; 2881 2882 PetscFunctionBegin; 2883 ierr = PetscUseMethod(A,"MatDenseRestoreColumn_C",(Mat,PetscScalar**),(A,vals));CHKERRQ(ierr); 2884 PetscFunctionReturn(0); 2885 } 2886