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