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