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