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