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_INTERN 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 PetscFunctionReturn(0); 721 } 722 723 static PetscErrorCode MatLUFactorSymbolic_SeqDense(Mat fact,Mat A,IS row,IS col,const MatFactorInfo *info) 724 { 725 PetscFunctionBegin; 726 fact->preallocated = PETSC_TRUE; 727 fact->assembled = PETSC_TRUE; 728 fact->ops->lufactornumeric = MatLUFactorNumeric_SeqDense; 729 PetscFunctionReturn(0); 730 } 731 732 PETSC_INTERN PetscErrorCode MatGetFactor_seqdense_petsc(Mat A,MatFactorType ftype,Mat *fact) 733 { 734 PetscErrorCode ierr; 735 736 PetscFunctionBegin; 737 ierr = MatCreate(PetscObjectComm((PetscObject)A),fact);CHKERRQ(ierr); 738 ierr = MatSetSizes(*fact,A->rmap->n,A->cmap->n,A->rmap->n,A->cmap->n);CHKERRQ(ierr); 739 ierr = MatSetType(*fact,((PetscObject)A)->type_name);CHKERRQ(ierr); 740 if (ftype == MAT_FACTOR_LU) { 741 (*fact)->ops->lufactorsymbolic = MatLUFactorSymbolic_SeqDense; 742 } else { 743 (*fact)->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_SeqDense; 744 } 745 (*fact)->factortype = ftype; 746 747 ierr = PetscFree((*fact)->solvertype);CHKERRQ(ierr); 748 ierr = PetscStrallocpy(MATSOLVERPETSC,&(*fact)->solvertype);CHKERRQ(ierr); 749 PetscFunctionReturn(0); 750 } 751 752 /* ------------------------------------------------------------------*/ 753 static PetscErrorCode MatSOR_SeqDense(Mat A,Vec bb,PetscReal omega,MatSORType flag,PetscReal shift,PetscInt its,PetscInt lits,Vec xx) 754 { 755 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 756 PetscScalar *x,*v = mat->v,zero = 0.0,xt; 757 const PetscScalar *b; 758 PetscErrorCode ierr; 759 PetscInt m = A->rmap->n,i; 760 PetscBLASInt o = 1,bm; 761 762 PetscFunctionBegin; 763 if (shift == -1) shift = 0.0; /* negative shift indicates do not error on zero diagonal; this code never zeros on zero diagonal */ 764 ierr = PetscBLASIntCast(m,&bm);CHKERRQ(ierr); 765 if (flag & SOR_ZERO_INITIAL_GUESS) { 766 /* this is a hack fix, should have another version without the second BLASdot */ 767 ierr = VecSet(xx,zero);CHKERRQ(ierr); 768 } 769 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 770 ierr = VecGetArrayRead(bb,&b);CHKERRQ(ierr); 771 its = its*lits; 772 if (its <= 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Relaxation requires global its %D and local its %D both positive",its,lits); 773 while (its--) { 774 if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP) { 775 for (i=0; i<m; i++) { 776 PetscStackCallBLAS("BLASdot",xt = b[i] - BLASdot_(&bm,v+i,&bm,x,&o)); 777 x[i] = (1. - omega)*x[i] + omega*(xt+v[i + i*m]*x[i])/(v[i + i*m]+shift); 778 } 779 } 780 if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP) { 781 for (i=m-1; i>=0; i--) { 782 PetscStackCallBLAS("BLASdot",xt = b[i] - BLASdot_(&bm,v+i,&bm,x,&o)); 783 x[i] = (1. - omega)*x[i] + omega*(xt+v[i + i*m]*x[i])/(v[i + i*m]+shift); 784 } 785 } 786 } 787 ierr = VecRestoreArrayRead(bb,&b);CHKERRQ(ierr); 788 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 789 PetscFunctionReturn(0); 790 } 791 792 /* -----------------------------------------------------------------*/ 793 PETSC_INTERN PetscErrorCode MatMultTranspose_SeqDense(Mat A,Vec xx,Vec yy) 794 { 795 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 796 const PetscScalar *v = mat->v,*x; 797 PetscScalar *y; 798 PetscErrorCode ierr; 799 PetscBLASInt m, n,_One=1; 800 PetscScalar _DOne=1.0,_DZero=0.0; 801 802 PetscFunctionBegin; 803 ierr = PetscBLASIntCast(A->rmap->n,&m);CHKERRQ(ierr); 804 ierr = PetscBLASIntCast(A->cmap->n,&n);CHKERRQ(ierr); 805 if (!A->rmap->n || !A->cmap->n) PetscFunctionReturn(0); 806 ierr = VecGetArrayRead(xx,&x);CHKERRQ(ierr); 807 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 808 PetscStackCallBLAS("BLASgemv",BLASgemv_("T",&m,&n,&_DOne,v,&mat->lda,x,&_One,&_DZero,y,&_One)); 809 ierr = VecRestoreArrayRead(xx,&x);CHKERRQ(ierr); 810 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 811 ierr = PetscLogFlops(2.0*A->rmap->n*A->cmap->n - A->cmap->n);CHKERRQ(ierr); 812 PetscFunctionReturn(0); 813 } 814 815 PETSC_INTERN PetscErrorCode MatMult_SeqDense(Mat A,Vec xx,Vec yy) 816 { 817 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 818 PetscScalar *y,_DOne=1.0,_DZero=0.0; 819 PetscErrorCode ierr; 820 PetscBLASInt m, n, _One=1; 821 const PetscScalar *v = mat->v,*x; 822 823 PetscFunctionBegin; 824 ierr = PetscBLASIntCast(A->rmap->n,&m);CHKERRQ(ierr); 825 ierr = PetscBLASIntCast(A->cmap->n,&n);CHKERRQ(ierr); 826 if (!A->rmap->n || !A->cmap->n) PetscFunctionReturn(0); 827 ierr = VecGetArrayRead(xx,&x);CHKERRQ(ierr); 828 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 829 PetscStackCallBLAS("BLASgemv",BLASgemv_("N",&m,&n,&_DOne,v,&(mat->lda),x,&_One,&_DZero,y,&_One)); 830 ierr = VecRestoreArrayRead(xx,&x);CHKERRQ(ierr); 831 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 832 ierr = PetscLogFlops(2.0*A->rmap->n*A->cmap->n - A->rmap->n);CHKERRQ(ierr); 833 PetscFunctionReturn(0); 834 } 835 836 PETSC_INTERN PetscErrorCode MatMultAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy) 837 { 838 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 839 const PetscScalar *v = mat->v,*x; 840 PetscScalar *y,_DOne=1.0; 841 PetscErrorCode ierr; 842 PetscBLASInt m, n, _One=1; 843 844 PetscFunctionBegin; 845 ierr = PetscBLASIntCast(A->rmap->n,&m);CHKERRQ(ierr); 846 ierr = PetscBLASIntCast(A->cmap->n,&n);CHKERRQ(ierr); 847 if (!A->rmap->n || !A->cmap->n) PetscFunctionReturn(0); 848 if (zz != yy) {ierr = VecCopy(zz,yy);CHKERRQ(ierr);} 849 ierr = VecGetArrayRead(xx,&x);CHKERRQ(ierr); 850 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 851 PetscStackCallBLAS("BLASgemv",BLASgemv_("N",&m,&n,&_DOne,v,&(mat->lda),x,&_One,&_DOne,y,&_One)); 852 ierr = VecRestoreArrayRead(xx,&x);CHKERRQ(ierr); 853 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 854 ierr = PetscLogFlops(2.0*A->rmap->n*A->cmap->n);CHKERRQ(ierr); 855 PetscFunctionReturn(0); 856 } 857 858 static PetscErrorCode MatMultTransposeAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy) 859 { 860 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 861 const PetscScalar *v = mat->v,*x; 862 PetscScalar *y; 863 PetscErrorCode ierr; 864 PetscBLASInt m, n, _One=1; 865 PetscScalar _DOne=1.0; 866 867 PetscFunctionBegin; 868 ierr = PetscBLASIntCast(A->rmap->n,&m);CHKERRQ(ierr); 869 ierr = PetscBLASIntCast(A->cmap->n,&n);CHKERRQ(ierr); 870 if (!A->rmap->n || !A->cmap->n) PetscFunctionReturn(0); 871 if (zz != yy) {ierr = VecCopy(zz,yy);CHKERRQ(ierr);} 872 ierr = VecGetArrayRead(xx,&x);CHKERRQ(ierr); 873 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 874 PetscStackCallBLAS("BLASgemv",BLASgemv_("T",&m,&n,&_DOne,v,&(mat->lda),x,&_One,&_DOne,y,&_One)); 875 ierr = VecRestoreArrayRead(xx,&x);CHKERRQ(ierr); 876 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 877 ierr = PetscLogFlops(2.0*A->rmap->n*A->cmap->n);CHKERRQ(ierr); 878 PetscFunctionReturn(0); 879 } 880 881 /* -----------------------------------------------------------------*/ 882 static PetscErrorCode MatGetRow_SeqDense(Mat A,PetscInt row,PetscInt *ncols,PetscInt **cols,PetscScalar **vals) 883 { 884 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 885 PetscScalar *v; 886 PetscErrorCode ierr; 887 PetscInt i; 888 889 PetscFunctionBegin; 890 *ncols = A->cmap->n; 891 if (cols) { 892 ierr = PetscMalloc1(A->cmap->n+1,cols);CHKERRQ(ierr); 893 for (i=0; i<A->cmap->n; i++) (*cols)[i] = i; 894 } 895 if (vals) { 896 ierr = PetscMalloc1(A->cmap->n+1,vals);CHKERRQ(ierr); 897 v = mat->v + row; 898 for (i=0; i<A->cmap->n; i++) {(*vals)[i] = *v; v += mat->lda;} 899 } 900 PetscFunctionReturn(0); 901 } 902 903 static PetscErrorCode MatRestoreRow_SeqDense(Mat A,PetscInt row,PetscInt *ncols,PetscInt **cols,PetscScalar **vals) 904 { 905 PetscErrorCode ierr; 906 907 PetscFunctionBegin; 908 if (cols) {ierr = PetscFree(*cols);CHKERRQ(ierr);} 909 if (vals) {ierr = PetscFree(*vals);CHKERRQ(ierr); } 910 PetscFunctionReturn(0); 911 } 912 /* ----------------------------------------------------------------*/ 913 static PetscErrorCode MatSetValues_SeqDense(Mat A,PetscInt m,const PetscInt indexm[],PetscInt n,const PetscInt indexn[],const PetscScalar v[],InsertMode addv) 914 { 915 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 916 PetscInt i,j,idx=0; 917 918 PetscFunctionBegin; 919 if (!mat->roworiented) { 920 if (addv == INSERT_VALUES) { 921 for (j=0; j<n; j++) { 922 if (indexn[j] < 0) {idx += m; continue;} 923 #if defined(PETSC_USE_DEBUG) 924 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); 925 #endif 926 for (i=0; i<m; i++) { 927 if (indexm[i] < 0) {idx++; continue;} 928 #if defined(PETSC_USE_DEBUG) 929 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); 930 #endif 931 mat->v[indexn[j]*mat->lda + indexm[i]] = v[idx++]; 932 } 933 } 934 } else { 935 for (j=0; j<n; j++) { 936 if (indexn[j] < 0) {idx += m; continue;} 937 #if defined(PETSC_USE_DEBUG) 938 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); 939 #endif 940 for (i=0; i<m; i++) { 941 if (indexm[i] < 0) {idx++; continue;} 942 #if defined(PETSC_USE_DEBUG) 943 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); 944 #endif 945 mat->v[indexn[j]*mat->lda + indexm[i]] += v[idx++]; 946 } 947 } 948 } 949 } else { 950 if (addv == INSERT_VALUES) { 951 for (i=0; i<m; i++) { 952 if (indexm[i] < 0) { idx += n; continue;} 953 #if defined(PETSC_USE_DEBUG) 954 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); 955 #endif 956 for (j=0; j<n; j++) { 957 if (indexn[j] < 0) { idx++; continue;} 958 #if defined(PETSC_USE_DEBUG) 959 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); 960 #endif 961 mat->v[indexn[j]*mat->lda + indexm[i]] = v[idx++]; 962 } 963 } 964 } else { 965 for (i=0; i<m; i++) { 966 if (indexm[i] < 0) { idx += n; continue;} 967 #if defined(PETSC_USE_DEBUG) 968 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); 969 #endif 970 for (j=0; j<n; j++) { 971 if (indexn[j] < 0) { idx++; continue;} 972 #if defined(PETSC_USE_DEBUG) 973 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); 974 #endif 975 mat->v[indexn[j]*mat->lda + indexm[i]] += v[idx++]; 976 } 977 } 978 } 979 } 980 PetscFunctionReturn(0); 981 } 982 983 static PetscErrorCode MatGetValues_SeqDense(Mat A,PetscInt m,const PetscInt indexm[],PetscInt n,const PetscInt indexn[],PetscScalar v[]) 984 { 985 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 986 PetscInt i,j; 987 988 PetscFunctionBegin; 989 /* row-oriented output */ 990 for (i=0; i<m; i++) { 991 if (indexm[i] < 0) {v += n;continue;} 992 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); 993 for (j=0; j<n; j++) { 994 if (indexn[j] < 0) {v++; continue;} 995 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); 996 *v++ = mat->v[indexn[j]*mat->lda + indexm[i]]; 997 } 998 } 999 PetscFunctionReturn(0); 1000 } 1001 1002 /* -----------------------------------------------------------------*/ 1003 1004 static PetscErrorCode MatLoad_SeqDense(Mat newmat,PetscViewer viewer) 1005 { 1006 Mat_SeqDense *a; 1007 PetscErrorCode ierr; 1008 PetscInt *scols,i,j,nz,header[4]; 1009 int fd; 1010 PetscMPIInt size; 1011 PetscInt *rowlengths = 0,M,N,*cols,grows,gcols; 1012 PetscScalar *vals,*svals,*v,*w; 1013 MPI_Comm comm; 1014 1015 PetscFunctionBegin; 1016 /* force binary viewer to load .info file if it has not yet done so */ 1017 ierr = PetscViewerSetUp(viewer);CHKERRQ(ierr); 1018 ierr = PetscObjectGetComm((PetscObject)viewer,&comm);CHKERRQ(ierr); 1019 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 1020 if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"view must have one processor"); 1021 ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 1022 ierr = PetscBinaryRead(fd,header,4,PETSC_INT);CHKERRQ(ierr); 1023 if (header[0] != MAT_FILE_CLASSID) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"Not matrix object"); 1024 M = header[1]; N = header[2]; nz = header[3]; 1025 1026 /* set global size if not set already*/ 1027 if (newmat->rmap->n < 0 && newmat->rmap->N < 0 && newmat->cmap->n < 0 && newmat->cmap->N < 0) { 1028 ierr = MatSetSizes(newmat,M,N,M,N);CHKERRQ(ierr); 1029 } else { 1030 /* if sizes and type are already set, check if the vector global sizes are correct */ 1031 ierr = MatGetSize(newmat,&grows,&gcols);CHKERRQ(ierr); 1032 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); 1033 } 1034 a = (Mat_SeqDense*)newmat->data; 1035 if (!a->user_alloc) { 1036 ierr = MatSeqDenseSetPreallocation(newmat,NULL);CHKERRQ(ierr); 1037 } 1038 1039 if (nz == MATRIX_BINARY_FORMAT_DENSE) { /* matrix in file is dense */ 1040 a = (Mat_SeqDense*)newmat->data; 1041 v = a->v; 1042 /* Allocate some temp space to read in the values and then flip them 1043 from row major to column major */ 1044 ierr = PetscMalloc1(M*N > 0 ? M*N : 1,&w);CHKERRQ(ierr); 1045 /* read in nonzero values */ 1046 ierr = PetscBinaryRead(fd,w,M*N,PETSC_SCALAR);CHKERRQ(ierr); 1047 /* now flip the values and store them in the matrix*/ 1048 for (j=0; j<N; j++) { 1049 for (i=0; i<M; i++) { 1050 *v++ =w[i*N+j]; 1051 } 1052 } 1053 ierr = PetscFree(w);CHKERRQ(ierr); 1054 } else { 1055 /* read row lengths */ 1056 ierr = PetscMalloc1(M+1,&rowlengths);CHKERRQ(ierr); 1057 ierr = PetscBinaryRead(fd,rowlengths,M,PETSC_INT);CHKERRQ(ierr); 1058 1059 a = (Mat_SeqDense*)newmat->data; 1060 v = a->v; 1061 1062 /* read column indices and nonzeros */ 1063 ierr = PetscMalloc1(nz+1,&scols);CHKERRQ(ierr); 1064 cols = scols; 1065 ierr = PetscBinaryRead(fd,cols,nz,PETSC_INT);CHKERRQ(ierr); 1066 ierr = PetscMalloc1(nz+1,&svals);CHKERRQ(ierr); 1067 vals = svals; 1068 ierr = PetscBinaryRead(fd,vals,nz,PETSC_SCALAR);CHKERRQ(ierr); 1069 1070 /* insert into matrix */ 1071 for (i=0; i<M; i++) { 1072 for (j=0; j<rowlengths[i]; j++) v[i+M*scols[j]] = svals[j]; 1073 svals += rowlengths[i]; scols += rowlengths[i]; 1074 } 1075 ierr = PetscFree(vals);CHKERRQ(ierr); 1076 ierr = PetscFree(cols);CHKERRQ(ierr); 1077 ierr = PetscFree(rowlengths);CHKERRQ(ierr); 1078 } 1079 ierr = MatAssemblyBegin(newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1080 ierr = MatAssemblyEnd(newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1081 PetscFunctionReturn(0); 1082 } 1083 1084 static PetscErrorCode MatView_SeqDense_ASCII(Mat A,PetscViewer viewer) 1085 { 1086 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1087 PetscErrorCode ierr; 1088 PetscInt i,j; 1089 const char *name; 1090 PetscScalar *v; 1091 PetscViewerFormat format; 1092 #if defined(PETSC_USE_COMPLEX) 1093 PetscBool allreal = PETSC_TRUE; 1094 #endif 1095 1096 PetscFunctionBegin; 1097 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 1098 if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) { 1099 PetscFunctionReturn(0); /* do nothing for now */ 1100 } else if (format == PETSC_VIEWER_ASCII_COMMON) { 1101 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);CHKERRQ(ierr); 1102 for (i=0; i<A->rmap->n; i++) { 1103 v = a->v + i; 1104 ierr = PetscViewerASCIIPrintf(viewer,"row %D:",i);CHKERRQ(ierr); 1105 for (j=0; j<A->cmap->n; j++) { 1106 #if defined(PETSC_USE_COMPLEX) 1107 if (PetscRealPart(*v) != 0.0 && PetscImaginaryPart(*v) != 0.0) { 1108 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g + %g i) ",j,(double)PetscRealPart(*v),(double)PetscImaginaryPart(*v));CHKERRQ(ierr); 1109 } else if (PetscRealPart(*v)) { 1110 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g) ",j,(double)PetscRealPart(*v));CHKERRQ(ierr); 1111 } 1112 #else 1113 if (*v) { 1114 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g) ",j,(double)*v);CHKERRQ(ierr); 1115 } 1116 #endif 1117 v += a->lda; 1118 } 1119 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 1120 } 1121 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);CHKERRQ(ierr); 1122 } else { 1123 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);CHKERRQ(ierr); 1124 #if defined(PETSC_USE_COMPLEX) 1125 /* determine if matrix has all real values */ 1126 v = a->v; 1127 for (i=0; i<A->rmap->n*A->cmap->n; i++) { 1128 if (PetscImaginaryPart(v[i])) { allreal = PETSC_FALSE; break;} 1129 } 1130 #endif 1131 if (format == PETSC_VIEWER_ASCII_MATLAB) { 1132 ierr = PetscObjectGetName((PetscObject)A,&name);CHKERRQ(ierr); 1133 ierr = PetscViewerASCIIPrintf(viewer,"%% Size = %D %D \n",A->rmap->n,A->cmap->n);CHKERRQ(ierr); 1134 ierr = PetscViewerASCIIPrintf(viewer,"%s = zeros(%D,%D);\n",name,A->rmap->n,A->cmap->n);CHKERRQ(ierr); 1135 ierr = PetscViewerASCIIPrintf(viewer,"%s = [\n",name);CHKERRQ(ierr); 1136 } 1137 1138 for (i=0; i<A->rmap->n; i++) { 1139 v = a->v + i; 1140 for (j=0; j<A->cmap->n; j++) { 1141 #if defined(PETSC_USE_COMPLEX) 1142 if (allreal) { 1143 ierr = PetscViewerASCIIPrintf(viewer,"%18.16e ",(double)PetscRealPart(*v));CHKERRQ(ierr); 1144 } else { 1145 ierr = PetscViewerASCIIPrintf(viewer,"%18.16e + %18.16ei ",(double)PetscRealPart(*v),(double)PetscImaginaryPart(*v));CHKERRQ(ierr); 1146 } 1147 #else 1148 ierr = PetscViewerASCIIPrintf(viewer,"%18.16e ",(double)*v);CHKERRQ(ierr); 1149 #endif 1150 v += a->lda; 1151 } 1152 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 1153 } 1154 if (format == PETSC_VIEWER_ASCII_MATLAB) { 1155 ierr = PetscViewerASCIIPrintf(viewer,"];\n");CHKERRQ(ierr); 1156 } 1157 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);CHKERRQ(ierr); 1158 } 1159 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 1160 PetscFunctionReturn(0); 1161 } 1162 1163 static PetscErrorCode MatView_SeqDense_Binary(Mat A,PetscViewer viewer) 1164 { 1165 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1166 PetscErrorCode ierr; 1167 int fd; 1168 PetscInt ict,j,n = A->cmap->n,m = A->rmap->n,i,*col_lens,nz = m*n; 1169 PetscScalar *v,*anonz,*vals; 1170 PetscViewerFormat format; 1171 1172 PetscFunctionBegin; 1173 ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 1174 1175 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 1176 if (format == PETSC_VIEWER_NATIVE) { 1177 /* store the matrix as a dense matrix */ 1178 ierr = PetscMalloc1(4,&col_lens);CHKERRQ(ierr); 1179 1180 col_lens[0] = MAT_FILE_CLASSID; 1181 col_lens[1] = m; 1182 col_lens[2] = n; 1183 col_lens[3] = MATRIX_BINARY_FORMAT_DENSE; 1184 1185 ierr = PetscBinaryWrite(fd,col_lens,4,PETSC_INT,PETSC_TRUE);CHKERRQ(ierr); 1186 ierr = PetscFree(col_lens);CHKERRQ(ierr); 1187 1188 /* write out matrix, by rows */ 1189 ierr = PetscMalloc1(m*n+1,&vals);CHKERRQ(ierr); 1190 v = a->v; 1191 for (j=0; j<n; j++) { 1192 for (i=0; i<m; i++) { 1193 vals[j + i*n] = *v++; 1194 } 1195 } 1196 ierr = PetscBinaryWrite(fd,vals,n*m,PETSC_SCALAR,PETSC_FALSE);CHKERRQ(ierr); 1197 ierr = PetscFree(vals);CHKERRQ(ierr); 1198 } else { 1199 ierr = PetscMalloc1(4+nz,&col_lens);CHKERRQ(ierr); 1200 1201 col_lens[0] = MAT_FILE_CLASSID; 1202 col_lens[1] = m; 1203 col_lens[2] = n; 1204 col_lens[3] = nz; 1205 1206 /* store lengths of each row and write (including header) to file */ 1207 for (i=0; i<m; i++) col_lens[4+i] = n; 1208 ierr = PetscBinaryWrite(fd,col_lens,4+m,PETSC_INT,PETSC_TRUE);CHKERRQ(ierr); 1209 1210 /* Possibly should write in smaller increments, not whole matrix at once? */ 1211 /* store column indices (zero start index) */ 1212 ict = 0; 1213 for (i=0; i<m; i++) { 1214 for (j=0; j<n; j++) col_lens[ict++] = j; 1215 } 1216 ierr = PetscBinaryWrite(fd,col_lens,nz,PETSC_INT,PETSC_FALSE);CHKERRQ(ierr); 1217 ierr = PetscFree(col_lens);CHKERRQ(ierr); 1218 1219 /* store nonzero values */ 1220 ierr = PetscMalloc1(nz+1,&anonz);CHKERRQ(ierr); 1221 ict = 0; 1222 for (i=0; i<m; i++) { 1223 v = a->v + i; 1224 for (j=0; j<n; j++) { 1225 anonz[ict++] = *v; v += a->lda; 1226 } 1227 } 1228 ierr = PetscBinaryWrite(fd,anonz,nz,PETSC_SCALAR,PETSC_FALSE);CHKERRQ(ierr); 1229 ierr = PetscFree(anonz);CHKERRQ(ierr); 1230 } 1231 PetscFunctionReturn(0); 1232 } 1233 1234 #include <petscdraw.h> 1235 static PetscErrorCode MatView_SeqDense_Draw_Zoom(PetscDraw draw,void *Aa) 1236 { 1237 Mat A = (Mat) Aa; 1238 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1239 PetscErrorCode ierr; 1240 PetscInt m = A->rmap->n,n = A->cmap->n,i,j; 1241 int color = PETSC_DRAW_WHITE; 1242 PetscScalar *v = a->v; 1243 PetscViewer viewer; 1244 PetscReal xl,yl,xr,yr,x_l,x_r,y_l,y_r; 1245 PetscViewerFormat format; 1246 1247 PetscFunctionBegin; 1248 ierr = PetscObjectQuery((PetscObject)A,"Zoomviewer",(PetscObject*)&viewer);CHKERRQ(ierr); 1249 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 1250 ierr = PetscDrawGetCoordinates(draw,&xl,&yl,&xr,&yr);CHKERRQ(ierr); 1251 1252 /* Loop over matrix elements drawing boxes */ 1253 1254 if (format != PETSC_VIEWER_DRAW_CONTOUR) { 1255 ierr = PetscDrawCollectiveBegin(draw);CHKERRQ(ierr); 1256 /* Blue for negative and Red for positive */ 1257 for (j = 0; j < n; j++) { 1258 x_l = j; x_r = x_l + 1.0; 1259 for (i = 0; i < m; i++) { 1260 y_l = m - i - 1.0; 1261 y_r = y_l + 1.0; 1262 if (PetscRealPart(v[j*m+i]) > 0.) { 1263 color = PETSC_DRAW_RED; 1264 } else if (PetscRealPart(v[j*m+i]) < 0.) { 1265 color = PETSC_DRAW_BLUE; 1266 } else { 1267 continue; 1268 } 1269 ierr = PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);CHKERRQ(ierr); 1270 } 1271 } 1272 ierr = PetscDrawCollectiveEnd(draw);CHKERRQ(ierr); 1273 } else { 1274 /* use contour shading to indicate magnitude of values */ 1275 /* first determine max of all nonzero values */ 1276 PetscReal minv = 0.0, maxv = 0.0; 1277 PetscDraw popup; 1278 1279 for (i=0; i < m*n; i++) { 1280 if (PetscAbsScalar(v[i]) > maxv) maxv = PetscAbsScalar(v[i]); 1281 } 1282 if (minv >= maxv) maxv = minv + PETSC_SMALL; 1283 ierr = PetscDrawGetPopup(draw,&popup);CHKERRQ(ierr); 1284 ierr = PetscDrawScalePopup(popup,minv,maxv);CHKERRQ(ierr); 1285 1286 ierr = PetscDrawCollectiveBegin(draw);CHKERRQ(ierr); 1287 for (j=0; j<n; j++) { 1288 x_l = j; 1289 x_r = x_l + 1.0; 1290 for (i=0; i<m; i++) { 1291 y_l = m - i - 1.0; 1292 y_r = y_l + 1.0; 1293 color = PetscDrawRealToColor(PetscAbsScalar(v[j*m+i]),minv,maxv); 1294 ierr = PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);CHKERRQ(ierr); 1295 } 1296 } 1297 ierr = PetscDrawCollectiveEnd(draw);CHKERRQ(ierr); 1298 } 1299 PetscFunctionReturn(0); 1300 } 1301 1302 static PetscErrorCode MatView_SeqDense_Draw(Mat A,PetscViewer viewer) 1303 { 1304 PetscDraw draw; 1305 PetscBool isnull; 1306 PetscReal xr,yr,xl,yl,h,w; 1307 PetscErrorCode ierr; 1308 1309 PetscFunctionBegin; 1310 ierr = PetscViewerDrawGetDraw(viewer,0,&draw);CHKERRQ(ierr); 1311 ierr = PetscDrawIsNull(draw,&isnull);CHKERRQ(ierr); 1312 if (isnull) PetscFunctionReturn(0); 1313 1314 xr = A->cmap->n; yr = A->rmap->n; h = yr/10.0; w = xr/10.0; 1315 xr += w; yr += h; xl = -w; yl = -h; 1316 ierr = PetscDrawSetCoordinates(draw,xl,yl,xr,yr);CHKERRQ(ierr); 1317 ierr = PetscObjectCompose((PetscObject)A,"Zoomviewer",(PetscObject)viewer);CHKERRQ(ierr); 1318 ierr = PetscDrawZoom(draw,MatView_SeqDense_Draw_Zoom,A);CHKERRQ(ierr); 1319 ierr = PetscObjectCompose((PetscObject)A,"Zoomviewer",NULL);CHKERRQ(ierr); 1320 ierr = PetscDrawSave(draw);CHKERRQ(ierr); 1321 PetscFunctionReturn(0); 1322 } 1323 1324 PetscErrorCode MatView_SeqDense(Mat A,PetscViewer viewer) 1325 { 1326 PetscErrorCode ierr; 1327 PetscBool iascii,isbinary,isdraw; 1328 1329 PetscFunctionBegin; 1330 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 1331 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);CHKERRQ(ierr); 1332 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);CHKERRQ(ierr); 1333 1334 if (iascii) { 1335 ierr = MatView_SeqDense_ASCII(A,viewer);CHKERRQ(ierr); 1336 } else if (isbinary) { 1337 ierr = MatView_SeqDense_Binary(A,viewer);CHKERRQ(ierr); 1338 } else if (isdraw) { 1339 ierr = MatView_SeqDense_Draw(A,viewer);CHKERRQ(ierr); 1340 } 1341 PetscFunctionReturn(0); 1342 } 1343 1344 static PetscErrorCode MatDestroy_SeqDense(Mat mat) 1345 { 1346 Mat_SeqDense *l = (Mat_SeqDense*)mat->data; 1347 PetscErrorCode ierr; 1348 1349 PetscFunctionBegin; 1350 #if defined(PETSC_USE_LOG) 1351 PetscLogObjectState((PetscObject)mat,"Rows %D Cols %D",mat->rmap->n,mat->cmap->n); 1352 #endif 1353 ierr = PetscFree(l->pivots);CHKERRQ(ierr); 1354 ierr = PetscFree(l->fwork);CHKERRQ(ierr); 1355 ierr = MatDestroy(&l->ptapwork);CHKERRQ(ierr); 1356 if (!l->user_alloc) {ierr = PetscFree(l->v);CHKERRQ(ierr);} 1357 ierr = PetscFree(mat->data);CHKERRQ(ierr); 1358 1359 ierr = PetscObjectChangeTypeName((PetscObject)mat,0);CHKERRQ(ierr); 1360 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetArray_C",NULL);CHKERRQ(ierr); 1361 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatDenseRestoreArray_C",NULL);CHKERRQ(ierr); 1362 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatConvert_seqdense_seqaij_C",NULL);CHKERRQ(ierr); 1363 #if defined(PETSC_HAVE_ELEMENTAL) 1364 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatConvert_seqdense_elemental_C",NULL);CHKERRQ(ierr); 1365 #endif 1366 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatSeqDenseSetPreallocation_C",NULL);CHKERRQ(ierr); 1367 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatMatMult_seqaij_seqdense_C",NULL);CHKERRQ(ierr); 1368 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatMatMultSymbolic_seqaij_seqdense_C",NULL);CHKERRQ(ierr); 1369 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatMatMultNumeric_seqaij_seqdense_C",NULL);CHKERRQ(ierr); 1370 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatPtAP_seqaij_seqdense_C",NULL);CHKERRQ(ierr); 1371 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatTransposeMatMult_seqaij_seqdense_C",NULL);CHKERRQ(ierr); 1372 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatTransposeMatMultSymbolic_seqaij_seqdense_C",NULL);CHKERRQ(ierr); 1373 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatTransposeMatMultNumeric_seqaij_seqdense_C",NULL);CHKERRQ(ierr); 1374 PetscFunctionReturn(0); 1375 } 1376 1377 static PetscErrorCode MatTranspose_SeqDense(Mat A,MatReuse reuse,Mat *matout) 1378 { 1379 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1380 PetscErrorCode ierr; 1381 PetscInt k,j,m,n,M; 1382 PetscScalar *v,tmp; 1383 1384 PetscFunctionBegin; 1385 v = mat->v; m = A->rmap->n; M = mat->lda; n = A->cmap->n; 1386 if (reuse == MAT_INPLACE_MATRIX) { /* in place transpose */ 1387 if (m != n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Can not transpose non-square matrix in place"); 1388 else { 1389 for (j=0; j<m; j++) { 1390 for (k=0; k<j; k++) { 1391 tmp = v[j + k*M]; 1392 v[j + k*M] = v[k + j*M]; 1393 v[k + j*M] = tmp; 1394 } 1395 } 1396 } 1397 } else { /* out-of-place transpose */ 1398 Mat tmat; 1399 Mat_SeqDense *tmatd; 1400 PetscScalar *v2; 1401 PetscInt M2; 1402 1403 if (reuse == MAT_INITIAL_MATRIX) { 1404 ierr = MatCreate(PetscObjectComm((PetscObject)A),&tmat);CHKERRQ(ierr); 1405 ierr = MatSetSizes(tmat,A->cmap->n,A->rmap->n,A->cmap->n,A->rmap->n);CHKERRQ(ierr); 1406 ierr = MatSetType(tmat,((PetscObject)A)->type_name);CHKERRQ(ierr); 1407 ierr = MatSeqDenseSetPreallocation(tmat,NULL);CHKERRQ(ierr); 1408 } else { 1409 tmat = *matout; 1410 } 1411 tmatd = (Mat_SeqDense*)tmat->data; 1412 v = mat->v; v2 = tmatd->v; M2 = tmatd->lda; 1413 for (j=0; j<n; j++) { 1414 for (k=0; k<m; k++) v2[j + k*M2] = v[k + j*M]; 1415 } 1416 ierr = MatAssemblyBegin(tmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1417 ierr = MatAssemblyEnd(tmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1418 *matout = tmat; 1419 } 1420 PetscFunctionReturn(0); 1421 } 1422 1423 static PetscErrorCode MatEqual_SeqDense(Mat A1,Mat A2,PetscBool *flg) 1424 { 1425 Mat_SeqDense *mat1 = (Mat_SeqDense*)A1->data; 1426 Mat_SeqDense *mat2 = (Mat_SeqDense*)A2->data; 1427 PetscInt i,j; 1428 PetscScalar *v1,*v2; 1429 1430 PetscFunctionBegin; 1431 if (A1->rmap->n != A2->rmap->n) {*flg = PETSC_FALSE; PetscFunctionReturn(0);} 1432 if (A1->cmap->n != A2->cmap->n) {*flg = PETSC_FALSE; PetscFunctionReturn(0);} 1433 for (i=0; i<A1->rmap->n; i++) { 1434 v1 = mat1->v+i; v2 = mat2->v+i; 1435 for (j=0; j<A1->cmap->n; j++) { 1436 if (*v1 != *v2) {*flg = PETSC_FALSE; PetscFunctionReturn(0);} 1437 v1 += mat1->lda; v2 += mat2->lda; 1438 } 1439 } 1440 *flg = PETSC_TRUE; 1441 PetscFunctionReturn(0); 1442 } 1443 1444 static PetscErrorCode MatGetDiagonal_SeqDense(Mat A,Vec v) 1445 { 1446 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1447 PetscErrorCode ierr; 1448 PetscInt i,n,len; 1449 PetscScalar *x,zero = 0.0; 1450 1451 PetscFunctionBegin; 1452 ierr = VecSet(v,zero);CHKERRQ(ierr); 1453 ierr = VecGetSize(v,&n);CHKERRQ(ierr); 1454 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 1455 len = PetscMin(A->rmap->n,A->cmap->n); 1456 if (n != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming mat and vec"); 1457 for (i=0; i<len; i++) { 1458 x[i] = mat->v[i*mat->lda + i]; 1459 } 1460 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 1461 PetscFunctionReturn(0); 1462 } 1463 1464 static PetscErrorCode MatDiagonalScale_SeqDense(Mat A,Vec ll,Vec rr) 1465 { 1466 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1467 const PetscScalar *l,*r; 1468 PetscScalar x,*v; 1469 PetscErrorCode ierr; 1470 PetscInt i,j,m = A->rmap->n,n = A->cmap->n; 1471 1472 PetscFunctionBegin; 1473 if (ll) { 1474 ierr = VecGetSize(ll,&m);CHKERRQ(ierr); 1475 ierr = VecGetArrayRead(ll,&l);CHKERRQ(ierr); 1476 if (m != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Left scaling vec wrong size"); 1477 for (i=0; i<m; i++) { 1478 x = l[i]; 1479 v = mat->v + i; 1480 for (j=0; j<n; j++) { (*v) *= x; v+= mat->lda;} 1481 } 1482 ierr = VecRestoreArrayRead(ll,&l);CHKERRQ(ierr); 1483 ierr = PetscLogFlops(1.0*n*m);CHKERRQ(ierr); 1484 } 1485 if (rr) { 1486 ierr = VecGetSize(rr,&n);CHKERRQ(ierr); 1487 ierr = VecGetArrayRead(rr,&r);CHKERRQ(ierr); 1488 if (n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Right scaling vec wrong size"); 1489 for (i=0; i<n; i++) { 1490 x = r[i]; 1491 v = mat->v + i*mat->lda; 1492 for (j=0; j<m; j++) (*v++) *= x; 1493 } 1494 ierr = VecRestoreArrayRead(rr,&r);CHKERRQ(ierr); 1495 ierr = PetscLogFlops(1.0*n*m);CHKERRQ(ierr); 1496 } 1497 PetscFunctionReturn(0); 1498 } 1499 1500 static PetscErrorCode MatNorm_SeqDense(Mat A,NormType type,PetscReal *nrm) 1501 { 1502 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1503 PetscScalar *v = mat->v; 1504 PetscReal sum = 0.0; 1505 PetscInt lda =mat->lda,m=A->rmap->n,i,j; 1506 PetscErrorCode ierr; 1507 1508 PetscFunctionBegin; 1509 if (type == NORM_FROBENIUS) { 1510 if (lda>m) { 1511 for (j=0; j<A->cmap->n; j++) { 1512 v = mat->v+j*lda; 1513 for (i=0; i<m; i++) { 1514 sum += PetscRealPart(PetscConj(*v)*(*v)); v++; 1515 } 1516 } 1517 } else { 1518 #if defined(PETSC_USE_REAL___FP16) 1519 PetscBLASInt one = 1,cnt = A->cmap->n*A->rmap->n; 1520 *nrm = BLASnrm2_(&cnt,v,&one); 1521 } 1522 #else 1523 for (i=0; i<A->cmap->n*A->rmap->n; i++) { 1524 sum += PetscRealPart(PetscConj(*v)*(*v)); v++; 1525 } 1526 } 1527 *nrm = PetscSqrtReal(sum); 1528 #endif 1529 ierr = PetscLogFlops(2.0*A->cmap->n*A->rmap->n);CHKERRQ(ierr); 1530 } else if (type == NORM_1) { 1531 *nrm = 0.0; 1532 for (j=0; j<A->cmap->n; j++) { 1533 v = mat->v + j*mat->lda; 1534 sum = 0.0; 1535 for (i=0; i<A->rmap->n; i++) { 1536 sum += PetscAbsScalar(*v); v++; 1537 } 1538 if (sum > *nrm) *nrm = sum; 1539 } 1540 ierr = PetscLogFlops(1.0*A->cmap->n*A->rmap->n);CHKERRQ(ierr); 1541 } else if (type == NORM_INFINITY) { 1542 *nrm = 0.0; 1543 for (j=0; j<A->rmap->n; j++) { 1544 v = mat->v + j; 1545 sum = 0.0; 1546 for (i=0; i<A->cmap->n; i++) { 1547 sum += PetscAbsScalar(*v); v += mat->lda; 1548 } 1549 if (sum > *nrm) *nrm = sum; 1550 } 1551 ierr = PetscLogFlops(1.0*A->cmap->n*A->rmap->n);CHKERRQ(ierr); 1552 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No two norm"); 1553 PetscFunctionReturn(0); 1554 } 1555 1556 static PetscErrorCode MatSetOption_SeqDense(Mat A,MatOption op,PetscBool flg) 1557 { 1558 Mat_SeqDense *aij = (Mat_SeqDense*)A->data; 1559 PetscErrorCode ierr; 1560 1561 PetscFunctionBegin; 1562 switch (op) { 1563 case MAT_ROW_ORIENTED: 1564 aij->roworiented = flg; 1565 break; 1566 case MAT_NEW_NONZERO_LOCATIONS: 1567 case MAT_NEW_NONZERO_LOCATION_ERR: 1568 case MAT_NEW_NONZERO_ALLOCATION_ERR: 1569 case MAT_NEW_DIAGONALS: 1570 case MAT_KEEP_NONZERO_PATTERN: 1571 case MAT_IGNORE_OFF_PROC_ENTRIES: 1572 case MAT_USE_HASH_TABLE: 1573 case MAT_IGNORE_LOWER_TRIANGULAR: 1574 ierr = PetscInfo1(A,"Option %s ignored\n",MatOptions[op]);CHKERRQ(ierr); 1575 break; 1576 case MAT_SPD: 1577 case MAT_SYMMETRIC: 1578 case MAT_STRUCTURALLY_SYMMETRIC: 1579 case MAT_HERMITIAN: 1580 case MAT_SYMMETRY_ETERNAL: 1581 /* These options are handled directly by MatSetOption() */ 1582 break; 1583 default: 1584 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"unknown option %s",MatOptions[op]); 1585 } 1586 PetscFunctionReturn(0); 1587 } 1588 1589 static PetscErrorCode MatZeroEntries_SeqDense(Mat A) 1590 { 1591 Mat_SeqDense *l = (Mat_SeqDense*)A->data; 1592 PetscErrorCode ierr; 1593 PetscInt lda=l->lda,m=A->rmap->n,j; 1594 1595 PetscFunctionBegin; 1596 if (lda>m) { 1597 for (j=0; j<A->cmap->n; j++) { 1598 ierr = PetscMemzero(l->v+j*lda,m*sizeof(PetscScalar));CHKERRQ(ierr); 1599 } 1600 } else { 1601 ierr = PetscMemzero(l->v,A->rmap->n*A->cmap->n*sizeof(PetscScalar));CHKERRQ(ierr); 1602 } 1603 PetscFunctionReturn(0); 1604 } 1605 1606 static PetscErrorCode MatZeroRows_SeqDense(Mat A,PetscInt N,const PetscInt rows[],PetscScalar diag,Vec x,Vec b) 1607 { 1608 PetscErrorCode ierr; 1609 Mat_SeqDense *l = (Mat_SeqDense*)A->data; 1610 PetscInt m = l->lda, n = A->cmap->n, i,j; 1611 PetscScalar *slot,*bb; 1612 const PetscScalar *xx; 1613 1614 PetscFunctionBegin; 1615 #if defined(PETSC_USE_DEBUG) 1616 for (i=0; i<N; i++) { 1617 if (rows[i] < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative row requested to be zeroed"); 1618 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); 1619 } 1620 #endif 1621 1622 /* fix right hand side if needed */ 1623 if (x && b) { 1624 ierr = VecGetArrayRead(x,&xx);CHKERRQ(ierr); 1625 ierr = VecGetArray(b,&bb);CHKERRQ(ierr); 1626 for (i=0; i<N; i++) bb[rows[i]] = diag*xx[rows[i]]; 1627 ierr = VecRestoreArrayRead(x,&xx);CHKERRQ(ierr); 1628 ierr = VecRestoreArray(b,&bb);CHKERRQ(ierr); 1629 } 1630 1631 for (i=0; i<N; i++) { 1632 slot = l->v + rows[i]; 1633 for (j=0; j<n; j++) { *slot = 0.0; slot += m;} 1634 } 1635 if (diag != 0.0) { 1636 if (A->rmap->n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only coded for square matrices"); 1637 for (i=0; i<N; i++) { 1638 slot = l->v + (m+1)*rows[i]; 1639 *slot = diag; 1640 } 1641 } 1642 PetscFunctionReturn(0); 1643 } 1644 1645 static PetscErrorCode MatDenseGetArray_SeqDense(Mat A,PetscScalar *array[]) 1646 { 1647 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1648 1649 PetscFunctionBegin; 1650 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"); 1651 *array = mat->v; 1652 PetscFunctionReturn(0); 1653 } 1654 1655 static PetscErrorCode MatDenseRestoreArray_SeqDense(Mat A,PetscScalar *array[]) 1656 { 1657 PetscFunctionBegin; 1658 *array = 0; /* user cannot accidently use the array later */ 1659 PetscFunctionReturn(0); 1660 } 1661 1662 /*@C 1663 MatDenseGetArray - gives access to the array where the data for a SeqDense matrix is stored 1664 1665 Not Collective 1666 1667 Input Parameter: 1668 . mat - a MATSEQDENSE or MATMPIDENSE matrix 1669 1670 Output Parameter: 1671 . array - pointer to the data 1672 1673 Level: intermediate 1674 1675 .seealso: MatDenseRestoreArray() 1676 @*/ 1677 PetscErrorCode MatDenseGetArray(Mat A,PetscScalar **array) 1678 { 1679 PetscErrorCode ierr; 1680 1681 PetscFunctionBegin; 1682 ierr = PetscUseMethod(A,"MatDenseGetArray_C",(Mat,PetscScalar**),(A,array));CHKERRQ(ierr); 1683 PetscFunctionReturn(0); 1684 } 1685 1686 /*@C 1687 MatDenseRestoreArray - returns access to the array where the data for a dense matrix is stored obtained by MatDenseGetArray() 1688 1689 Not Collective 1690 1691 Input Parameters: 1692 . mat - a MATSEQDENSE or MATMPIDENSE matrix 1693 . array - pointer to the data 1694 1695 Level: intermediate 1696 1697 .seealso: MatDenseGetArray() 1698 @*/ 1699 PetscErrorCode MatDenseRestoreArray(Mat A,PetscScalar **array) 1700 { 1701 PetscErrorCode ierr; 1702 1703 PetscFunctionBegin; 1704 ierr = PetscUseMethod(A,"MatDenseRestoreArray_C",(Mat,PetscScalar**),(A,array));CHKERRQ(ierr); 1705 PetscFunctionReturn(0); 1706 } 1707 1708 static PetscErrorCode MatCreateSubMatrix_SeqDense(Mat A,IS isrow,IS iscol,PetscInt cs,MatReuse scall,Mat *B) 1709 { 1710 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1711 PetscErrorCode ierr; 1712 PetscInt i,j,nrows,ncols; 1713 const PetscInt *irow,*icol; 1714 PetscScalar *av,*bv,*v = mat->v; 1715 Mat newmat; 1716 1717 PetscFunctionBegin; 1718 ierr = ISGetIndices(isrow,&irow);CHKERRQ(ierr); 1719 ierr = ISGetIndices(iscol,&icol);CHKERRQ(ierr); 1720 ierr = ISGetLocalSize(isrow,&nrows);CHKERRQ(ierr); 1721 ierr = ISGetLocalSize(iscol,&ncols);CHKERRQ(ierr); 1722 1723 /* Check submatrixcall */ 1724 if (scall == MAT_REUSE_MATRIX) { 1725 PetscInt n_cols,n_rows; 1726 ierr = MatGetSize(*B,&n_rows,&n_cols);CHKERRQ(ierr); 1727 if (n_rows != nrows || n_cols != ncols) { 1728 /* resize the result matrix to match number of requested rows/columns */ 1729 ierr = MatSetSizes(*B,nrows,ncols,nrows,ncols);CHKERRQ(ierr); 1730 } 1731 newmat = *B; 1732 } else { 1733 /* Create and fill new matrix */ 1734 ierr = MatCreate(PetscObjectComm((PetscObject)A),&newmat);CHKERRQ(ierr); 1735 ierr = MatSetSizes(newmat,nrows,ncols,nrows,ncols);CHKERRQ(ierr); 1736 ierr = MatSetType(newmat,((PetscObject)A)->type_name);CHKERRQ(ierr); 1737 ierr = MatSeqDenseSetPreallocation(newmat,NULL);CHKERRQ(ierr); 1738 } 1739 1740 /* Now extract the data pointers and do the copy,column at a time */ 1741 bv = ((Mat_SeqDense*)newmat->data)->v; 1742 1743 for (i=0; i<ncols; i++) { 1744 av = v + mat->lda*icol[i]; 1745 for (j=0; j<nrows; j++) *bv++ = av[irow[j]]; 1746 } 1747 1748 /* Assemble the matrices so that the correct flags are set */ 1749 ierr = MatAssemblyBegin(newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1750 ierr = MatAssemblyEnd(newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1751 1752 /* Free work space */ 1753 ierr = ISRestoreIndices(isrow,&irow);CHKERRQ(ierr); 1754 ierr = ISRestoreIndices(iscol,&icol);CHKERRQ(ierr); 1755 *B = newmat; 1756 PetscFunctionReturn(0); 1757 } 1758 1759 static PetscErrorCode MatCreateSubMatrices_SeqDense(Mat A,PetscInt n,const IS irow[],const IS icol[],MatReuse scall,Mat *B[]) 1760 { 1761 PetscErrorCode ierr; 1762 PetscInt i; 1763 1764 PetscFunctionBegin; 1765 if (scall == MAT_INITIAL_MATRIX) { 1766 ierr = PetscCalloc1(n+1,B);CHKERRQ(ierr); 1767 } 1768 1769 for (i=0; i<n; i++) { 1770 ierr = MatCreateSubMatrix_SeqDense(A,irow[i],icol[i],PETSC_DECIDE,scall,&(*B)[i]);CHKERRQ(ierr); 1771 } 1772 PetscFunctionReturn(0); 1773 } 1774 1775 static PetscErrorCode MatAssemblyBegin_SeqDense(Mat mat,MatAssemblyType mode) 1776 { 1777 PetscFunctionBegin; 1778 PetscFunctionReturn(0); 1779 } 1780 1781 static PetscErrorCode MatAssemblyEnd_SeqDense(Mat mat,MatAssemblyType mode) 1782 { 1783 PetscFunctionBegin; 1784 PetscFunctionReturn(0); 1785 } 1786 1787 static PetscErrorCode MatCopy_SeqDense(Mat A,Mat B,MatStructure str) 1788 { 1789 Mat_SeqDense *a = (Mat_SeqDense*)A->data,*b = (Mat_SeqDense*)B->data; 1790 PetscErrorCode ierr; 1791 PetscInt lda1=a->lda,lda2=b->lda, m=A->rmap->n,n=A->cmap->n, j; 1792 1793 PetscFunctionBegin; 1794 /* If the two matrices don't have the same copy implementation, they aren't compatible for fast copy. */ 1795 if (A->ops->copy != B->ops->copy) { 1796 ierr = MatCopy_Basic(A,B,str);CHKERRQ(ierr); 1797 PetscFunctionReturn(0); 1798 } 1799 if (m != B->rmap->n || n != B->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"size(B) != size(A)"); 1800 if (lda1>m || lda2>m) { 1801 for (j=0; j<n; j++) { 1802 ierr = PetscMemcpy(b->v+j*lda2,a->v+j*lda1,m*sizeof(PetscScalar));CHKERRQ(ierr); 1803 } 1804 } else { 1805 ierr = PetscMemcpy(b->v,a->v,A->rmap->n*A->cmap->n*sizeof(PetscScalar));CHKERRQ(ierr); 1806 } 1807 PetscFunctionReturn(0); 1808 } 1809 1810 static PetscErrorCode MatSetUp_SeqDense(Mat A) 1811 { 1812 PetscErrorCode ierr; 1813 1814 PetscFunctionBegin; 1815 ierr = MatSeqDenseSetPreallocation(A,0);CHKERRQ(ierr); 1816 PetscFunctionReturn(0); 1817 } 1818 1819 static PetscErrorCode MatConjugate_SeqDense(Mat A) 1820 { 1821 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1822 PetscInt i,nz = A->rmap->n*A->cmap->n; 1823 PetscScalar *aa = a->v; 1824 1825 PetscFunctionBegin; 1826 for (i=0; i<nz; i++) aa[i] = PetscConj(aa[i]); 1827 PetscFunctionReturn(0); 1828 } 1829 1830 static PetscErrorCode MatRealPart_SeqDense(Mat A) 1831 { 1832 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1833 PetscInt i,nz = A->rmap->n*A->cmap->n; 1834 PetscScalar *aa = a->v; 1835 1836 PetscFunctionBegin; 1837 for (i=0; i<nz; i++) aa[i] = PetscRealPart(aa[i]); 1838 PetscFunctionReturn(0); 1839 } 1840 1841 static PetscErrorCode MatImaginaryPart_SeqDense(Mat A) 1842 { 1843 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1844 PetscInt i,nz = A->rmap->n*A->cmap->n; 1845 PetscScalar *aa = a->v; 1846 1847 PetscFunctionBegin; 1848 for (i=0; i<nz; i++) aa[i] = PetscImaginaryPart(aa[i]); 1849 PetscFunctionReturn(0); 1850 } 1851 1852 /* ----------------------------------------------------------------*/ 1853 PETSC_INTERN PetscErrorCode MatMatMult_SeqDense_SeqDense(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C) 1854 { 1855 PetscErrorCode ierr; 1856 1857 PetscFunctionBegin; 1858 if (scall == MAT_INITIAL_MATRIX) { 1859 ierr = PetscLogEventBegin(MAT_MatMultSymbolic,A,B,0,0);CHKERRQ(ierr); 1860 ierr = MatMatMultSymbolic_SeqDense_SeqDense(A,B,fill,C);CHKERRQ(ierr); 1861 ierr = PetscLogEventEnd(MAT_MatMultSymbolic,A,B,0,0);CHKERRQ(ierr); 1862 } 1863 ierr = PetscLogEventBegin(MAT_MatMultNumeric,A,B,0,0);CHKERRQ(ierr); 1864 ierr = MatMatMultNumeric_SeqDense_SeqDense(A,B,*C);CHKERRQ(ierr); 1865 ierr = PetscLogEventEnd(MAT_MatMultNumeric,A,B,0,0);CHKERRQ(ierr); 1866 PetscFunctionReturn(0); 1867 } 1868 1869 PetscErrorCode MatMatMultSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat *C) 1870 { 1871 PetscErrorCode ierr; 1872 PetscInt m=A->rmap->n,n=B->cmap->n; 1873 Mat Cmat; 1874 1875 PetscFunctionBegin; 1876 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); 1877 ierr = MatCreate(PETSC_COMM_SELF,&Cmat);CHKERRQ(ierr); 1878 ierr = MatSetSizes(Cmat,m,n,m,n);CHKERRQ(ierr); 1879 ierr = MatSetType(Cmat,MATSEQDENSE);CHKERRQ(ierr); 1880 ierr = MatSeqDenseSetPreallocation(Cmat,NULL);CHKERRQ(ierr); 1881 1882 *C = Cmat; 1883 PetscFunctionReturn(0); 1884 } 1885 1886 PetscErrorCode MatMatMultNumeric_SeqDense_SeqDense(Mat A,Mat B,Mat C) 1887 { 1888 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1889 Mat_SeqDense *b = (Mat_SeqDense*)B->data; 1890 Mat_SeqDense *c = (Mat_SeqDense*)C->data; 1891 PetscBLASInt m,n,k; 1892 PetscScalar _DOne=1.0,_DZero=0.0; 1893 PetscErrorCode ierr; 1894 PetscBool flg; 1895 1896 PetscFunctionBegin; 1897 ierr = PetscObjectTypeCompare((PetscObject)B,MATSEQDENSE,&flg);CHKERRQ(ierr); 1898 if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Second matrix must be dense"); 1899 1900 /* Handle case where where user provided the final C matrix rather than calling MatMatMult() with MAT_INITIAL_MATRIX*/ 1901 ierr = PetscObjectTypeCompare((PetscObject)A,MATSEQAIJ,&flg);CHKERRQ(ierr); 1902 if (flg) { 1903 C->ops->matmultnumeric = MatMatMultNumeric_SeqAIJ_SeqDense; 1904 ierr = (*C->ops->matmultnumeric)(A,B,C);CHKERRQ(ierr); 1905 PetscFunctionReturn(0); 1906 } 1907 1908 ierr = PetscObjectTypeCompare((PetscObject)A,MATSEQDENSE,&flg);CHKERRQ(ierr); 1909 if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"First matrix must be dense"); 1910 ierr = PetscBLASIntCast(C->rmap->n,&m);CHKERRQ(ierr); 1911 ierr = PetscBLASIntCast(C->cmap->n,&n);CHKERRQ(ierr); 1912 ierr = PetscBLASIntCast(A->cmap->n,&k);CHKERRQ(ierr); 1913 PetscStackCallBLAS("BLASgemm",BLASgemm_("N","N",&m,&n,&k,&_DOne,a->v,&a->lda,b->v,&b->lda,&_DZero,c->v,&c->lda)); 1914 PetscFunctionReturn(0); 1915 } 1916 1917 PetscErrorCode MatTransposeMatMult_SeqDense_SeqDense(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C) 1918 { 1919 PetscErrorCode ierr; 1920 1921 PetscFunctionBegin; 1922 if (scall == MAT_INITIAL_MATRIX) { 1923 ierr = PetscLogEventBegin(MAT_TransposeMatMultSymbolic,A,B,0,0);CHKERRQ(ierr); 1924 ierr = MatTransposeMatMultSymbolic_SeqDense_SeqDense(A,B,fill,C);CHKERRQ(ierr); 1925 ierr = PetscLogEventEnd(MAT_TransposeMatMultSymbolic,A,B,0,0);CHKERRQ(ierr); 1926 } 1927 ierr = PetscLogEventBegin(MAT_TransposeMatMultNumeric,A,B,0,0);CHKERRQ(ierr); 1928 ierr = MatTransposeMatMultNumeric_SeqDense_SeqDense(A,B,*C);CHKERRQ(ierr); 1929 ierr = PetscLogEventEnd(MAT_TransposeMatMultNumeric,A,B,0,0);CHKERRQ(ierr); 1930 PetscFunctionReturn(0); 1931 } 1932 1933 PetscErrorCode MatTransposeMatMultSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat *C) 1934 { 1935 PetscErrorCode ierr; 1936 PetscInt m=A->cmap->n,n=B->cmap->n; 1937 Mat Cmat; 1938 1939 PetscFunctionBegin; 1940 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); 1941 ierr = MatCreate(PETSC_COMM_SELF,&Cmat);CHKERRQ(ierr); 1942 ierr = MatSetSizes(Cmat,m,n,m,n);CHKERRQ(ierr); 1943 ierr = MatSetType(Cmat,MATSEQDENSE);CHKERRQ(ierr); 1944 ierr = MatSeqDenseSetPreallocation(Cmat,NULL);CHKERRQ(ierr); 1945 1946 Cmat->assembled = PETSC_TRUE; 1947 1948 *C = Cmat; 1949 PetscFunctionReturn(0); 1950 } 1951 1952 PetscErrorCode MatTransposeMatMultNumeric_SeqDense_SeqDense(Mat A,Mat B,Mat C) 1953 { 1954 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1955 Mat_SeqDense *b = (Mat_SeqDense*)B->data; 1956 Mat_SeqDense *c = (Mat_SeqDense*)C->data; 1957 PetscBLASInt m,n,k; 1958 PetscScalar _DOne=1.0,_DZero=0.0; 1959 PetscErrorCode ierr; 1960 1961 PetscFunctionBegin; 1962 ierr = PetscBLASIntCast(C->rmap->n,&m);CHKERRQ(ierr); 1963 ierr = PetscBLASIntCast(C->cmap->n,&n);CHKERRQ(ierr); 1964 ierr = PetscBLASIntCast(A->rmap->n,&k);CHKERRQ(ierr); 1965 PetscStackCallBLAS("BLASgemm",BLASgemm_("T","N",&m,&n,&k,&_DOne,a->v,&a->lda,b->v,&b->lda,&_DZero,c->v,&c->lda)); 1966 PetscFunctionReturn(0); 1967 } 1968 1969 static PetscErrorCode MatGetRowMax_SeqDense(Mat A,Vec v,PetscInt idx[]) 1970 { 1971 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1972 PetscErrorCode ierr; 1973 PetscInt i,j,m = A->rmap->n,n = A->cmap->n,p; 1974 PetscScalar *x; 1975 MatScalar *aa = a->v; 1976 1977 PetscFunctionBegin; 1978 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 1979 1980 ierr = VecSet(v,0.0);CHKERRQ(ierr); 1981 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 1982 ierr = VecGetLocalSize(v,&p);CHKERRQ(ierr); 1983 if (p != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 1984 for (i=0; i<m; i++) { 1985 x[i] = aa[i]; if (idx) idx[i] = 0; 1986 for (j=1; j<n; j++) { 1987 if (PetscRealPart(x[i]) < PetscRealPart(aa[i+m*j])) {x[i] = aa[i + m*j]; if (idx) idx[i] = j;} 1988 } 1989 } 1990 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 1991 PetscFunctionReturn(0); 1992 } 1993 1994 static PetscErrorCode MatGetRowMaxAbs_SeqDense(Mat A,Vec v,PetscInt idx[]) 1995 { 1996 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1997 PetscErrorCode ierr; 1998 PetscInt i,j,m = A->rmap->n,n = A->cmap->n,p; 1999 PetscScalar *x; 2000 PetscReal atmp; 2001 MatScalar *aa = a->v; 2002 2003 PetscFunctionBegin; 2004 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2005 2006 ierr = VecSet(v,0.0);CHKERRQ(ierr); 2007 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2008 ierr = VecGetLocalSize(v,&p);CHKERRQ(ierr); 2009 if (p != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 2010 for (i=0; i<m; i++) { 2011 x[i] = PetscAbsScalar(aa[i]); 2012 for (j=1; j<n; j++) { 2013 atmp = PetscAbsScalar(aa[i+m*j]); 2014 if (PetscAbsScalar(x[i]) < atmp) {x[i] = atmp; if (idx) idx[i] = j;} 2015 } 2016 } 2017 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2018 PetscFunctionReturn(0); 2019 } 2020 2021 static PetscErrorCode MatGetRowMin_SeqDense(Mat A,Vec v,PetscInt idx[]) 2022 { 2023 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 2024 PetscErrorCode ierr; 2025 PetscInt i,j,m = A->rmap->n,n = A->cmap->n,p; 2026 PetscScalar *x; 2027 MatScalar *aa = a->v; 2028 2029 PetscFunctionBegin; 2030 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2031 2032 ierr = VecSet(v,0.0);CHKERRQ(ierr); 2033 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2034 ierr = VecGetLocalSize(v,&p);CHKERRQ(ierr); 2035 if (p != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 2036 for (i=0; i<m; i++) { 2037 x[i] = aa[i]; if (idx) idx[i] = 0; 2038 for (j=1; j<n; j++) { 2039 if (PetscRealPart(x[i]) > PetscRealPart(aa[i+m*j])) {x[i] = aa[i + m*j]; if (idx) idx[i] = j;} 2040 } 2041 } 2042 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2043 PetscFunctionReturn(0); 2044 } 2045 2046 static PetscErrorCode MatGetColumnVector_SeqDense(Mat A,Vec v,PetscInt col) 2047 { 2048 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 2049 PetscErrorCode ierr; 2050 PetscScalar *x; 2051 2052 PetscFunctionBegin; 2053 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2054 2055 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2056 ierr = PetscMemcpy(x,a->v+col*a->lda,A->rmap->n*sizeof(PetscScalar));CHKERRQ(ierr); 2057 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2058 PetscFunctionReturn(0); 2059 } 2060 2061 2062 PetscErrorCode MatGetColumnNorms_SeqDense(Mat A,NormType type,PetscReal *norms) 2063 { 2064 PetscErrorCode ierr; 2065 PetscInt i,j,m,n; 2066 PetscScalar *a; 2067 2068 PetscFunctionBegin; 2069 ierr = MatGetSize(A,&m,&n);CHKERRQ(ierr); 2070 ierr = PetscMemzero(norms,n*sizeof(PetscReal));CHKERRQ(ierr); 2071 ierr = MatDenseGetArray(A,&a);CHKERRQ(ierr); 2072 if (type == NORM_2) { 2073 for (i=0; i<n; i++) { 2074 for (j=0; j<m; j++) { 2075 norms[i] += PetscAbsScalar(a[j]*a[j]); 2076 } 2077 a += m; 2078 } 2079 } else if (type == NORM_1) { 2080 for (i=0; i<n; i++) { 2081 for (j=0; j<m; j++) { 2082 norms[i] += PetscAbsScalar(a[j]); 2083 } 2084 a += m; 2085 } 2086 } else if (type == NORM_INFINITY) { 2087 for (i=0; i<n; i++) { 2088 for (j=0; j<m; j++) { 2089 norms[i] = PetscMax(PetscAbsScalar(a[j]),norms[i]); 2090 } 2091 a += m; 2092 } 2093 } else SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Unknown NormType"); 2094 ierr = MatDenseRestoreArray(A,&a);CHKERRQ(ierr); 2095 if (type == NORM_2) { 2096 for (i=0; i<n; i++) norms[i] = PetscSqrtReal(norms[i]); 2097 } 2098 PetscFunctionReturn(0); 2099 } 2100 2101 static PetscErrorCode MatSetRandom_SeqDense(Mat x,PetscRandom rctx) 2102 { 2103 PetscErrorCode ierr; 2104 PetscScalar *a; 2105 PetscInt m,n,i; 2106 2107 PetscFunctionBegin; 2108 ierr = MatGetSize(x,&m,&n);CHKERRQ(ierr); 2109 ierr = MatDenseGetArray(x,&a);CHKERRQ(ierr); 2110 for (i=0; i<m*n; i++) { 2111 ierr = PetscRandomGetValue(rctx,a+i);CHKERRQ(ierr); 2112 } 2113 ierr = MatDenseRestoreArray(x,&a);CHKERRQ(ierr); 2114 PetscFunctionReturn(0); 2115 } 2116 2117 static PetscErrorCode MatMissingDiagonal_SeqDense(Mat A,PetscBool *missing,PetscInt *d) 2118 { 2119 PetscFunctionBegin; 2120 *missing = PETSC_FALSE; 2121 PetscFunctionReturn(0); 2122 } 2123 2124 2125 /* -------------------------------------------------------------------*/ 2126 static struct _MatOps MatOps_Values = { MatSetValues_SeqDense, 2127 MatGetRow_SeqDense, 2128 MatRestoreRow_SeqDense, 2129 MatMult_SeqDense, 2130 /* 4*/ MatMultAdd_SeqDense, 2131 MatMultTranspose_SeqDense, 2132 MatMultTransposeAdd_SeqDense, 2133 0, 2134 0, 2135 0, 2136 /* 10*/ 0, 2137 MatLUFactor_SeqDense, 2138 MatCholeskyFactor_SeqDense, 2139 MatSOR_SeqDense, 2140 MatTranspose_SeqDense, 2141 /* 15*/ MatGetInfo_SeqDense, 2142 MatEqual_SeqDense, 2143 MatGetDiagonal_SeqDense, 2144 MatDiagonalScale_SeqDense, 2145 MatNorm_SeqDense, 2146 /* 20*/ MatAssemblyBegin_SeqDense, 2147 MatAssemblyEnd_SeqDense, 2148 MatSetOption_SeqDense, 2149 MatZeroEntries_SeqDense, 2150 /* 24*/ MatZeroRows_SeqDense, 2151 0, 2152 0, 2153 0, 2154 0, 2155 /* 29*/ MatSetUp_SeqDense, 2156 0, 2157 0, 2158 0, 2159 0, 2160 /* 34*/ MatDuplicate_SeqDense, 2161 0, 2162 0, 2163 0, 2164 0, 2165 /* 39*/ MatAXPY_SeqDense, 2166 MatCreateSubMatrices_SeqDense, 2167 0, 2168 MatGetValues_SeqDense, 2169 MatCopy_SeqDense, 2170 /* 44*/ MatGetRowMax_SeqDense, 2171 MatScale_SeqDense, 2172 MatShift_Basic, 2173 0, 2174 MatZeroRowsColumns_SeqDense, 2175 /* 49*/ MatSetRandom_SeqDense, 2176 0, 2177 0, 2178 0, 2179 0, 2180 /* 54*/ 0, 2181 0, 2182 0, 2183 0, 2184 0, 2185 /* 59*/ 0, 2186 MatDestroy_SeqDense, 2187 MatView_SeqDense, 2188 0, 2189 0, 2190 /* 64*/ 0, 2191 0, 2192 0, 2193 0, 2194 0, 2195 /* 69*/ MatGetRowMaxAbs_SeqDense, 2196 0, 2197 0, 2198 0, 2199 0, 2200 /* 74*/ 0, 2201 0, 2202 0, 2203 0, 2204 0, 2205 /* 79*/ 0, 2206 0, 2207 0, 2208 0, 2209 /* 83*/ MatLoad_SeqDense, 2210 0, 2211 MatIsHermitian_SeqDense, 2212 0, 2213 0, 2214 0, 2215 /* 89*/ MatMatMult_SeqDense_SeqDense, 2216 MatMatMultSymbolic_SeqDense_SeqDense, 2217 MatMatMultNumeric_SeqDense_SeqDense, 2218 MatPtAP_SeqDense_SeqDense, 2219 MatPtAPSymbolic_SeqDense_SeqDense, 2220 /* 94*/ MatPtAPNumeric_SeqDense_SeqDense, 2221 0, 2222 0, 2223 0, 2224 0, 2225 /* 99*/ 0, 2226 0, 2227 0, 2228 MatConjugate_SeqDense, 2229 0, 2230 /*104*/ 0, 2231 MatRealPart_SeqDense, 2232 MatImaginaryPart_SeqDense, 2233 0, 2234 0, 2235 /*109*/ 0, 2236 0, 2237 MatGetRowMin_SeqDense, 2238 MatGetColumnVector_SeqDense, 2239 MatMissingDiagonal_SeqDense, 2240 /*114*/ 0, 2241 0, 2242 0, 2243 0, 2244 0, 2245 /*119*/ 0, 2246 0, 2247 0, 2248 0, 2249 0, 2250 /*124*/ 0, 2251 MatGetColumnNorms_SeqDense, 2252 0, 2253 0, 2254 0, 2255 /*129*/ 0, 2256 MatTransposeMatMult_SeqDense_SeqDense, 2257 MatTransposeMatMultSymbolic_SeqDense_SeqDense, 2258 MatTransposeMatMultNumeric_SeqDense_SeqDense, 2259 0, 2260 /*134*/ 0, 2261 0, 2262 0, 2263 0, 2264 0, 2265 /*139*/ 0, 2266 0, 2267 0 2268 }; 2269 2270 /*@C 2271 MatCreateSeqDense - Creates a sequential dense matrix that 2272 is stored in column major order (the usual Fortran 77 manner). Many 2273 of the matrix operations use the BLAS and LAPACK routines. 2274 2275 Collective on MPI_Comm 2276 2277 Input Parameters: 2278 + comm - MPI communicator, set to PETSC_COMM_SELF 2279 . m - number of rows 2280 . n - number of columns 2281 - data - optional location of matrix data in column major order. Set data=NULL for PETSc 2282 to control all matrix memory allocation. 2283 2284 Output Parameter: 2285 . A - the matrix 2286 2287 Notes: 2288 The data input variable is intended primarily for Fortran programmers 2289 who wish to allocate their own matrix memory space. Most users should 2290 set data=NULL. 2291 2292 Level: intermediate 2293 2294 .keywords: dense, matrix, LAPACK, BLAS 2295 2296 .seealso: MatCreate(), MatCreateDense(), MatSetValues() 2297 @*/ 2298 PetscErrorCode MatCreateSeqDense(MPI_Comm comm,PetscInt m,PetscInt n,PetscScalar *data,Mat *A) 2299 { 2300 PetscErrorCode ierr; 2301 2302 PetscFunctionBegin; 2303 ierr = MatCreate(comm,A);CHKERRQ(ierr); 2304 ierr = MatSetSizes(*A,m,n,m,n);CHKERRQ(ierr); 2305 ierr = MatSetType(*A,MATSEQDENSE);CHKERRQ(ierr); 2306 ierr = MatSeqDenseSetPreallocation(*A,data);CHKERRQ(ierr); 2307 PetscFunctionReturn(0); 2308 } 2309 2310 /*@C 2311 MatSeqDenseSetPreallocation - Sets the array used for storing the matrix elements 2312 2313 Collective on MPI_Comm 2314 2315 Input Parameters: 2316 + B - the matrix 2317 - data - the array (or NULL) 2318 2319 Notes: 2320 The data input variable is intended primarily for Fortran programmers 2321 who wish to allocate their own matrix memory space. Most users should 2322 need not call this routine. 2323 2324 Level: intermediate 2325 2326 .keywords: dense, matrix, LAPACK, BLAS 2327 2328 .seealso: MatCreate(), MatCreateDense(), MatSetValues(), MatSeqDenseSetLDA() 2329 2330 @*/ 2331 PetscErrorCode MatSeqDenseSetPreallocation(Mat B,PetscScalar data[]) 2332 { 2333 PetscErrorCode ierr; 2334 2335 PetscFunctionBegin; 2336 ierr = PetscTryMethod(B,"MatSeqDenseSetPreallocation_C",(Mat,PetscScalar[]),(B,data));CHKERRQ(ierr); 2337 PetscFunctionReturn(0); 2338 } 2339 2340 PetscErrorCode MatSeqDenseSetPreallocation_SeqDense(Mat B,PetscScalar *data) 2341 { 2342 Mat_SeqDense *b; 2343 PetscErrorCode ierr; 2344 2345 PetscFunctionBegin; 2346 B->preallocated = PETSC_TRUE; 2347 2348 ierr = PetscLayoutSetUp(B->rmap);CHKERRQ(ierr); 2349 ierr = PetscLayoutSetUp(B->cmap);CHKERRQ(ierr); 2350 2351 b = (Mat_SeqDense*)B->data; 2352 b->Mmax = B->rmap->n; 2353 b->Nmax = B->cmap->n; 2354 if (b->lda <= 0 || b->changelda) b->lda = B->rmap->n; 2355 2356 ierr = PetscIntMultError(b->lda,b->Nmax,NULL);CHKERRQ(ierr); 2357 if (!data) { /* petsc-allocated storage */ 2358 if (!b->user_alloc) { ierr = PetscFree(b->v);CHKERRQ(ierr); } 2359 ierr = PetscCalloc1((size_t)b->lda*b->Nmax,&b->v);CHKERRQ(ierr); 2360 ierr = PetscLogObjectMemory((PetscObject)B,b->lda*b->Nmax*sizeof(PetscScalar));CHKERRQ(ierr); 2361 2362 b->user_alloc = PETSC_FALSE; 2363 } else { /* user-allocated storage */ 2364 if (!b->user_alloc) { ierr = PetscFree(b->v);CHKERRQ(ierr); } 2365 b->v = data; 2366 b->user_alloc = PETSC_TRUE; 2367 } 2368 B->assembled = PETSC_TRUE; 2369 PetscFunctionReturn(0); 2370 } 2371 2372 #if defined(PETSC_HAVE_ELEMENTAL) 2373 PETSC_INTERN PetscErrorCode MatConvert_SeqDense_Elemental(Mat A, MatType newtype,MatReuse reuse,Mat *newmat) 2374 { 2375 Mat mat_elemental; 2376 PetscErrorCode ierr; 2377 PetscScalar *array,*v_colwise; 2378 PetscInt M=A->rmap->N,N=A->cmap->N,i,j,k,*rows,*cols; 2379 2380 PetscFunctionBegin; 2381 ierr = PetscMalloc3(M*N,&v_colwise,M,&rows,N,&cols);CHKERRQ(ierr); 2382 ierr = MatDenseGetArray(A,&array);CHKERRQ(ierr); 2383 /* convert column-wise array into row-wise v_colwise, see MatSetValues_Elemental() */ 2384 k = 0; 2385 for (j=0; j<N; j++) { 2386 cols[j] = j; 2387 for (i=0; i<M; i++) { 2388 v_colwise[j*M+i] = array[k++]; 2389 } 2390 } 2391 for (i=0; i<M; i++) { 2392 rows[i] = i; 2393 } 2394 ierr = MatDenseRestoreArray(A,&array);CHKERRQ(ierr); 2395 2396 ierr = MatCreate(PetscObjectComm((PetscObject)A), &mat_elemental);CHKERRQ(ierr); 2397 ierr = MatSetSizes(mat_elemental,PETSC_DECIDE,PETSC_DECIDE,M,N);CHKERRQ(ierr); 2398 ierr = MatSetType(mat_elemental,MATELEMENTAL);CHKERRQ(ierr); 2399 ierr = MatSetUp(mat_elemental);CHKERRQ(ierr); 2400 2401 /* PETSc-Elemental interaface uses axpy for setting off-processor entries, only ADD_VALUES is allowed */ 2402 ierr = MatSetValues(mat_elemental,M,rows,N,cols,v_colwise,ADD_VALUES);CHKERRQ(ierr); 2403 ierr = MatAssemblyBegin(mat_elemental, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2404 ierr = MatAssemblyEnd(mat_elemental, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2405 ierr = PetscFree3(v_colwise,rows,cols);CHKERRQ(ierr); 2406 2407 if (reuse == MAT_INPLACE_MATRIX) { 2408 ierr = MatHeaderReplace(A,&mat_elemental);CHKERRQ(ierr); 2409 } else { 2410 *newmat = mat_elemental; 2411 } 2412 PetscFunctionReturn(0); 2413 } 2414 #endif 2415 2416 /*@C 2417 MatSeqDenseSetLDA - Declare the leading dimension of the user-provided array 2418 2419 Input parameter: 2420 + A - the matrix 2421 - lda - the leading dimension 2422 2423 Notes: 2424 This routine is to be used in conjunction with MatSeqDenseSetPreallocation(); 2425 it asserts that the preallocation has a leading dimension (the LDA parameter 2426 of Blas and Lapack fame) larger than M, the first dimension of the matrix. 2427 2428 Level: intermediate 2429 2430 .keywords: dense, matrix, LAPACK, BLAS 2431 2432 .seealso: MatCreate(), MatCreateSeqDense(), MatSeqDenseSetPreallocation(), MatSetMaximumSize() 2433 2434 @*/ 2435 PetscErrorCode MatSeqDenseSetLDA(Mat B,PetscInt lda) 2436 { 2437 Mat_SeqDense *b = (Mat_SeqDense*)B->data; 2438 2439 PetscFunctionBegin; 2440 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); 2441 b->lda = lda; 2442 b->changelda = PETSC_FALSE; 2443 b->Mmax = PetscMax(b->Mmax,lda); 2444 PetscFunctionReturn(0); 2445 } 2446 2447 /*MC 2448 MATSEQDENSE - MATSEQDENSE = "seqdense" - A matrix type to be used for sequential dense matrices. 2449 2450 Options Database Keys: 2451 . -mat_type seqdense - sets the matrix type to "seqdense" during a call to MatSetFromOptions() 2452 2453 Level: beginner 2454 2455 .seealso: MatCreateSeqDense() 2456 2457 M*/ 2458 2459 PETSC_EXTERN PetscErrorCode MatCreate_SeqDense(Mat B) 2460 { 2461 Mat_SeqDense *b; 2462 PetscErrorCode ierr; 2463 PetscMPIInt size; 2464 2465 PetscFunctionBegin; 2466 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)B),&size);CHKERRQ(ierr); 2467 if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Comm must be of size 1"); 2468 2469 ierr = PetscNewLog(B,&b);CHKERRQ(ierr); 2470 ierr = PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 2471 B->data = (void*)b; 2472 2473 b->roworiented = PETSC_TRUE; 2474 2475 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseGetArray_C",MatDenseGetArray_SeqDense);CHKERRQ(ierr); 2476 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreArray_C",MatDenseRestoreArray_SeqDense);CHKERRQ(ierr); 2477 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqdense_seqaij_C",MatConvert_SeqDense_SeqAIJ);CHKERRQ(ierr); 2478 #if defined(PETSC_HAVE_ELEMENTAL) 2479 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqdense_elemental_C",MatConvert_SeqDense_Elemental);CHKERRQ(ierr); 2480 #endif 2481 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqDenseSetPreallocation_C",MatSeqDenseSetPreallocation_SeqDense);CHKERRQ(ierr); 2482 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMult_seqaij_seqdense_C",MatMatMult_SeqAIJ_SeqDense);CHKERRQ(ierr); 2483 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMultSymbolic_seqaij_seqdense_C",MatMatMultSymbolic_SeqAIJ_SeqDense);CHKERRQ(ierr); 2484 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMultNumeric_seqaij_seqdense_C",MatMatMultNumeric_SeqAIJ_SeqDense);CHKERRQ(ierr); 2485 ierr = PetscObjectComposeFunction((PetscObject)B,"MatPtAP_seqaij_seqdense_C",MatPtAP_SeqDense_SeqDense);CHKERRQ(ierr); 2486 ierr = PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMult_seqaij_seqdense_C",MatTransposeMatMult_SeqAIJ_SeqDense);CHKERRQ(ierr); 2487 ierr = PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultSymbolic_seqaij_seqdense_C",MatTransposeMatMultSymbolic_SeqAIJ_SeqDense);CHKERRQ(ierr); 2488 ierr = PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultNumeric_seqaij_seqdense_C",MatTransposeMatMultNumeric_SeqAIJ_SeqDense);CHKERRQ(ierr); 2489 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQDENSE);CHKERRQ(ierr); 2490 PetscFunctionReturn(0); 2491 } 2492