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