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 MatDensePlaceArray_SeqDense(Mat A,const PetscScalar array[]) 1240 { 1241 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1242 1243 PetscFunctionBegin; 1244 a->unplacedarray = a->v; 1245 a->unplaced_user_alloc = a->user_alloc; 1246 a->v = (PetscScalar*) array; 1247 PetscFunctionReturn(0); 1248 } 1249 1250 static PetscErrorCode MatDenseResetArray_SeqDense(Mat A) 1251 { 1252 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1253 1254 PetscFunctionBegin; 1255 a->v = a->unplacedarray; 1256 a->user_alloc = a->unplaced_user_alloc; 1257 a->unplacedarray = NULL; 1258 PetscFunctionReturn(0); 1259 } 1260 1261 static PetscErrorCode MatDestroy_SeqDense(Mat mat) 1262 { 1263 Mat_SeqDense *l = (Mat_SeqDense*)mat->data; 1264 PetscErrorCode ierr; 1265 1266 PetscFunctionBegin; 1267 #if defined(PETSC_USE_LOG) 1268 PetscLogObjectState((PetscObject)mat,"Rows %D Cols %D",mat->rmap->n,mat->cmap->n); 1269 #endif 1270 ierr = PetscFree(l->pivots);CHKERRQ(ierr); 1271 ierr = MatDestroy(&l->ptapwork);CHKERRQ(ierr); 1272 if (!l->user_alloc) {ierr = PetscFree(l->v);CHKERRQ(ierr);} 1273 ierr = PetscFree(mat->data);CHKERRQ(ierr); 1274 1275 ierr = PetscObjectChangeTypeName((PetscObject)mat,0);CHKERRQ(ierr); 1276 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetArray_C",NULL);CHKERRQ(ierr); 1277 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatDensePlaceArray_C",NULL);CHKERRQ(ierr); 1278 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatDenseResetArray_C",NULL);CHKERRQ(ierr); 1279 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatDenseRestoreArray_C",NULL);CHKERRQ(ierr); 1280 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatConvert_seqdense_seqaij_C",NULL);CHKERRQ(ierr); 1281 #if defined(PETSC_HAVE_ELEMENTAL) 1282 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatConvert_seqdense_elemental_C",NULL);CHKERRQ(ierr); 1283 #endif 1284 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatSeqDenseSetPreallocation_C",NULL);CHKERRQ(ierr); 1285 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatMatMult_seqaij_seqdense_C",NULL);CHKERRQ(ierr); 1286 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatMatMultSymbolic_seqaij_seqdense_C",NULL);CHKERRQ(ierr); 1287 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatMatMultNumeric_seqaij_seqdense_C",NULL);CHKERRQ(ierr); 1288 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatPtAP_seqaij_seqdense_C",NULL);CHKERRQ(ierr); 1289 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatTransposeMatMult_seqaij_seqdense_C",NULL);CHKERRQ(ierr); 1290 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatTransposeMatMultSymbolic_seqaij_seqdense_C",NULL);CHKERRQ(ierr); 1291 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatTransposeMatMultNumeric_seqaij_seqdense_C",NULL);CHKERRQ(ierr); 1292 PetscFunctionReturn(0); 1293 } 1294 1295 static PetscErrorCode MatTranspose_SeqDense(Mat A,MatReuse reuse,Mat *matout) 1296 { 1297 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1298 PetscErrorCode ierr; 1299 PetscInt k,j,m,n,M; 1300 PetscScalar *v,tmp; 1301 1302 PetscFunctionBegin; 1303 v = mat->v; m = A->rmap->n; M = mat->lda; n = A->cmap->n; 1304 if (reuse == MAT_INPLACE_MATRIX) { /* in place transpose */ 1305 if (m != n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Can not transpose non-square matrix in place"); 1306 else { 1307 for (j=0; j<m; j++) { 1308 for (k=0; k<j; k++) { 1309 tmp = v[j + k*M]; 1310 v[j + k*M] = v[k + j*M]; 1311 v[k + j*M] = tmp; 1312 } 1313 } 1314 } 1315 } else { /* out-of-place transpose */ 1316 Mat tmat; 1317 Mat_SeqDense *tmatd; 1318 PetscScalar *v2; 1319 PetscInt M2; 1320 1321 if (reuse == MAT_INITIAL_MATRIX) { 1322 ierr = MatCreate(PetscObjectComm((PetscObject)A),&tmat);CHKERRQ(ierr); 1323 ierr = MatSetSizes(tmat,A->cmap->n,A->rmap->n,A->cmap->n,A->rmap->n);CHKERRQ(ierr); 1324 ierr = MatSetType(tmat,((PetscObject)A)->type_name);CHKERRQ(ierr); 1325 ierr = MatSeqDenseSetPreallocation(tmat,NULL);CHKERRQ(ierr); 1326 } else { 1327 tmat = *matout; 1328 } 1329 tmatd = (Mat_SeqDense*)tmat->data; 1330 v = mat->v; v2 = tmatd->v; M2 = tmatd->lda; 1331 for (j=0; j<n; j++) { 1332 for (k=0; k<m; k++) v2[j + k*M2] = v[k + j*M]; 1333 } 1334 ierr = MatAssemblyBegin(tmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1335 ierr = MatAssemblyEnd(tmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1336 *matout = tmat; 1337 } 1338 PetscFunctionReturn(0); 1339 } 1340 1341 static PetscErrorCode MatEqual_SeqDense(Mat A1,Mat A2,PetscBool *flg) 1342 { 1343 Mat_SeqDense *mat1 = (Mat_SeqDense*)A1->data; 1344 Mat_SeqDense *mat2 = (Mat_SeqDense*)A2->data; 1345 PetscInt i,j; 1346 PetscScalar *v1,*v2; 1347 1348 PetscFunctionBegin; 1349 if (A1->rmap->n != A2->rmap->n) {*flg = PETSC_FALSE; PetscFunctionReturn(0);} 1350 if (A1->cmap->n != A2->cmap->n) {*flg = PETSC_FALSE; PetscFunctionReturn(0);} 1351 for (i=0; i<A1->rmap->n; i++) { 1352 v1 = mat1->v+i; v2 = mat2->v+i; 1353 for (j=0; j<A1->cmap->n; j++) { 1354 if (*v1 != *v2) {*flg = PETSC_FALSE; PetscFunctionReturn(0);} 1355 v1 += mat1->lda; v2 += mat2->lda; 1356 } 1357 } 1358 *flg = PETSC_TRUE; 1359 PetscFunctionReturn(0); 1360 } 1361 1362 static PetscErrorCode MatGetDiagonal_SeqDense(Mat A,Vec v) 1363 { 1364 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1365 PetscErrorCode ierr; 1366 PetscInt i,n,len; 1367 PetscScalar *x,zero = 0.0; 1368 1369 PetscFunctionBegin; 1370 ierr = VecSet(v,zero);CHKERRQ(ierr); 1371 ierr = VecGetSize(v,&n);CHKERRQ(ierr); 1372 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 1373 len = PetscMin(A->rmap->n,A->cmap->n); 1374 if (n != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming mat and vec"); 1375 for (i=0; i<len; i++) { 1376 x[i] = mat->v[i*mat->lda + i]; 1377 } 1378 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 1379 PetscFunctionReturn(0); 1380 } 1381 1382 static PetscErrorCode MatDiagonalScale_SeqDense(Mat A,Vec ll,Vec rr) 1383 { 1384 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1385 const PetscScalar *l,*r; 1386 PetscScalar x,*v; 1387 PetscErrorCode ierr; 1388 PetscInt i,j,m = A->rmap->n,n = A->cmap->n; 1389 1390 PetscFunctionBegin; 1391 if (ll) { 1392 ierr = VecGetSize(ll,&m);CHKERRQ(ierr); 1393 ierr = VecGetArrayRead(ll,&l);CHKERRQ(ierr); 1394 if (m != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Left scaling vec wrong size"); 1395 for (i=0; i<m; i++) { 1396 x = l[i]; 1397 v = mat->v + i; 1398 for (j=0; j<n; j++) { (*v) *= x; v+= mat->lda;} 1399 } 1400 ierr = VecRestoreArrayRead(ll,&l);CHKERRQ(ierr); 1401 ierr = PetscLogFlops(1.0*n*m);CHKERRQ(ierr); 1402 } 1403 if (rr) { 1404 ierr = VecGetSize(rr,&n);CHKERRQ(ierr); 1405 ierr = VecGetArrayRead(rr,&r);CHKERRQ(ierr); 1406 if (n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Right scaling vec wrong size"); 1407 for (i=0; i<n; i++) { 1408 x = r[i]; 1409 v = mat->v + i*mat->lda; 1410 for (j=0; j<m; j++) (*v++) *= x; 1411 } 1412 ierr = VecRestoreArrayRead(rr,&r);CHKERRQ(ierr); 1413 ierr = PetscLogFlops(1.0*n*m);CHKERRQ(ierr); 1414 } 1415 PetscFunctionReturn(0); 1416 } 1417 1418 static PetscErrorCode MatNorm_SeqDense(Mat A,NormType type,PetscReal *nrm) 1419 { 1420 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1421 PetscScalar *v = mat->v; 1422 PetscReal sum = 0.0; 1423 PetscInt lda =mat->lda,m=A->rmap->n,i,j; 1424 PetscErrorCode ierr; 1425 1426 PetscFunctionBegin; 1427 if (type == NORM_FROBENIUS) { 1428 if (lda>m) { 1429 for (j=0; j<A->cmap->n; j++) { 1430 v = mat->v+j*lda; 1431 for (i=0; i<m; i++) { 1432 sum += PetscRealPart(PetscConj(*v)*(*v)); v++; 1433 } 1434 } 1435 } else { 1436 #if defined(PETSC_USE_REAL___FP16) 1437 PetscBLASInt one = 1,cnt = A->cmap->n*A->rmap->n; 1438 *nrm = BLASnrm2_(&cnt,v,&one); 1439 } 1440 #else 1441 for (i=0; i<A->cmap->n*A->rmap->n; i++) { 1442 sum += PetscRealPart(PetscConj(*v)*(*v)); v++; 1443 } 1444 } 1445 *nrm = PetscSqrtReal(sum); 1446 #endif 1447 ierr = PetscLogFlops(2.0*A->cmap->n*A->rmap->n);CHKERRQ(ierr); 1448 } else if (type == NORM_1) { 1449 *nrm = 0.0; 1450 for (j=0; j<A->cmap->n; j++) { 1451 v = mat->v + j*mat->lda; 1452 sum = 0.0; 1453 for (i=0; i<A->rmap->n; i++) { 1454 sum += PetscAbsScalar(*v); v++; 1455 } 1456 if (sum > *nrm) *nrm = sum; 1457 } 1458 ierr = PetscLogFlops(1.0*A->cmap->n*A->rmap->n);CHKERRQ(ierr); 1459 } else if (type == NORM_INFINITY) { 1460 *nrm = 0.0; 1461 for (j=0; j<A->rmap->n; j++) { 1462 v = mat->v + j; 1463 sum = 0.0; 1464 for (i=0; i<A->cmap->n; i++) { 1465 sum += PetscAbsScalar(*v); v += mat->lda; 1466 } 1467 if (sum > *nrm) *nrm = sum; 1468 } 1469 ierr = PetscLogFlops(1.0*A->cmap->n*A->rmap->n);CHKERRQ(ierr); 1470 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No two norm"); 1471 PetscFunctionReturn(0); 1472 } 1473 1474 static PetscErrorCode MatSetOption_SeqDense(Mat A,MatOption op,PetscBool flg) 1475 { 1476 Mat_SeqDense *aij = (Mat_SeqDense*)A->data; 1477 PetscErrorCode ierr; 1478 1479 PetscFunctionBegin; 1480 switch (op) { 1481 case MAT_ROW_ORIENTED: 1482 aij->roworiented = flg; 1483 break; 1484 case MAT_NEW_NONZERO_LOCATIONS: 1485 case MAT_NEW_NONZERO_LOCATION_ERR: 1486 case MAT_NEW_NONZERO_ALLOCATION_ERR: 1487 case MAT_NEW_DIAGONALS: 1488 case MAT_KEEP_NONZERO_PATTERN: 1489 case MAT_IGNORE_OFF_PROC_ENTRIES: 1490 case MAT_USE_HASH_TABLE: 1491 case MAT_IGNORE_LOWER_TRIANGULAR: 1492 ierr = PetscInfo1(A,"Option %s ignored\n",MatOptions[op]);CHKERRQ(ierr); 1493 break; 1494 case MAT_SPD: 1495 case MAT_SYMMETRIC: 1496 case MAT_STRUCTURALLY_SYMMETRIC: 1497 case MAT_HERMITIAN: 1498 case MAT_SYMMETRY_ETERNAL: 1499 /* These options are handled directly by MatSetOption() */ 1500 break; 1501 default: 1502 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"unknown option %s",MatOptions[op]); 1503 } 1504 PetscFunctionReturn(0); 1505 } 1506 1507 static PetscErrorCode MatZeroEntries_SeqDense(Mat A) 1508 { 1509 Mat_SeqDense *l = (Mat_SeqDense*)A->data; 1510 PetscErrorCode ierr; 1511 PetscInt lda=l->lda,m=A->rmap->n,j; 1512 1513 PetscFunctionBegin; 1514 if (lda>m) { 1515 for (j=0; j<A->cmap->n; j++) { 1516 ierr = PetscMemzero(l->v+j*lda,m*sizeof(PetscScalar));CHKERRQ(ierr); 1517 } 1518 } else { 1519 ierr = PetscMemzero(l->v,A->rmap->n*A->cmap->n*sizeof(PetscScalar));CHKERRQ(ierr); 1520 } 1521 PetscFunctionReturn(0); 1522 } 1523 1524 static PetscErrorCode MatZeroRows_SeqDense(Mat A,PetscInt N,const PetscInt rows[],PetscScalar diag,Vec x,Vec b) 1525 { 1526 PetscErrorCode ierr; 1527 Mat_SeqDense *l = (Mat_SeqDense*)A->data; 1528 PetscInt m = l->lda, n = A->cmap->n, i,j; 1529 PetscScalar *slot,*bb; 1530 const PetscScalar *xx; 1531 1532 PetscFunctionBegin; 1533 #if defined(PETSC_USE_DEBUG) 1534 for (i=0; i<N; i++) { 1535 if (rows[i] < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative row requested to be zeroed"); 1536 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); 1537 } 1538 #endif 1539 1540 /* fix right hand side if needed */ 1541 if (x && b) { 1542 ierr = VecGetArrayRead(x,&xx);CHKERRQ(ierr); 1543 ierr = VecGetArray(b,&bb);CHKERRQ(ierr); 1544 for (i=0; i<N; i++) bb[rows[i]] = diag*xx[rows[i]]; 1545 ierr = VecRestoreArrayRead(x,&xx);CHKERRQ(ierr); 1546 ierr = VecRestoreArray(b,&bb);CHKERRQ(ierr); 1547 } 1548 1549 for (i=0; i<N; i++) { 1550 slot = l->v + rows[i]; 1551 for (j=0; j<n; j++) { *slot = 0.0; slot += m;} 1552 } 1553 if (diag != 0.0) { 1554 if (A->rmap->n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only coded for square matrices"); 1555 for (i=0; i<N; i++) { 1556 slot = l->v + (m+1)*rows[i]; 1557 *slot = diag; 1558 } 1559 } 1560 PetscFunctionReturn(0); 1561 } 1562 1563 static PetscErrorCode MatDenseGetArray_SeqDense(Mat A,PetscScalar *array[]) 1564 { 1565 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1566 1567 PetscFunctionBegin; 1568 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"); 1569 *array = mat->v; 1570 PetscFunctionReturn(0); 1571 } 1572 1573 static PetscErrorCode MatDenseRestoreArray_SeqDense(Mat A,PetscScalar *array[]) 1574 { 1575 PetscFunctionBegin; 1576 *array = 0; /* user cannot accidently use the array later */ 1577 PetscFunctionReturn(0); 1578 } 1579 1580 /*@C 1581 MatDenseGetArray - gives access to the array where the data for a SeqDense matrix is stored 1582 1583 Not Collective 1584 1585 Input Parameter: 1586 . mat - a MATSEQDENSE or MATMPIDENSE matrix 1587 1588 Output Parameter: 1589 . array - pointer to the data 1590 1591 Level: intermediate 1592 1593 .seealso: MatDenseRestoreArray() 1594 @*/ 1595 PetscErrorCode MatDenseGetArray(Mat A,PetscScalar **array) 1596 { 1597 PetscErrorCode ierr; 1598 1599 PetscFunctionBegin; 1600 ierr = PetscUseMethod(A,"MatDenseGetArray_C",(Mat,PetscScalar**),(A,array));CHKERRQ(ierr); 1601 PetscFunctionReturn(0); 1602 } 1603 1604 /*@C 1605 MatDenseRestoreArray - returns access to the array where the data for a dense matrix is stored obtained by MatDenseGetArray() 1606 1607 Not Collective 1608 1609 Input Parameters: 1610 . mat - a MATSEQDENSE or MATMPIDENSE matrix 1611 . array - pointer to the data 1612 1613 Level: intermediate 1614 1615 .seealso: MatDenseGetArray() 1616 @*/ 1617 PetscErrorCode MatDenseRestoreArray(Mat A,PetscScalar **array) 1618 { 1619 PetscErrorCode ierr; 1620 1621 PetscFunctionBegin; 1622 ierr = PetscUseMethod(A,"MatDenseRestoreArray_C",(Mat,PetscScalar**),(A,array));CHKERRQ(ierr); 1623 PetscFunctionReturn(0); 1624 } 1625 1626 static PetscErrorCode MatCreateSubMatrix_SeqDense(Mat A,IS isrow,IS iscol,PetscInt cs,MatReuse scall,Mat *B) 1627 { 1628 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1629 PetscErrorCode ierr; 1630 PetscInt i,j,nrows,ncols; 1631 const PetscInt *irow,*icol; 1632 PetscScalar *av,*bv,*v = mat->v; 1633 Mat newmat; 1634 1635 PetscFunctionBegin; 1636 ierr = ISGetIndices(isrow,&irow);CHKERRQ(ierr); 1637 ierr = ISGetIndices(iscol,&icol);CHKERRQ(ierr); 1638 ierr = ISGetLocalSize(isrow,&nrows);CHKERRQ(ierr); 1639 ierr = ISGetLocalSize(iscol,&ncols);CHKERRQ(ierr); 1640 1641 /* Check submatrixcall */ 1642 if (scall == MAT_REUSE_MATRIX) { 1643 PetscInt n_cols,n_rows; 1644 ierr = MatGetSize(*B,&n_rows,&n_cols);CHKERRQ(ierr); 1645 if (n_rows != nrows || n_cols != ncols) { 1646 /* resize the result matrix to match number of requested rows/columns */ 1647 ierr = MatSetSizes(*B,nrows,ncols,nrows,ncols);CHKERRQ(ierr); 1648 } 1649 newmat = *B; 1650 } else { 1651 /* Create and fill new matrix */ 1652 ierr = MatCreate(PetscObjectComm((PetscObject)A),&newmat);CHKERRQ(ierr); 1653 ierr = MatSetSizes(newmat,nrows,ncols,nrows,ncols);CHKERRQ(ierr); 1654 ierr = MatSetType(newmat,((PetscObject)A)->type_name);CHKERRQ(ierr); 1655 ierr = MatSeqDenseSetPreallocation(newmat,NULL);CHKERRQ(ierr); 1656 } 1657 1658 /* Now extract the data pointers and do the copy,column at a time */ 1659 bv = ((Mat_SeqDense*)newmat->data)->v; 1660 1661 for (i=0; i<ncols; i++) { 1662 av = v + mat->lda*icol[i]; 1663 for (j=0; j<nrows; j++) *bv++ = av[irow[j]]; 1664 } 1665 1666 /* Assemble the matrices so that the correct flags are set */ 1667 ierr = MatAssemblyBegin(newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1668 ierr = MatAssemblyEnd(newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1669 1670 /* Free work space */ 1671 ierr = ISRestoreIndices(isrow,&irow);CHKERRQ(ierr); 1672 ierr = ISRestoreIndices(iscol,&icol);CHKERRQ(ierr); 1673 *B = newmat; 1674 PetscFunctionReturn(0); 1675 } 1676 1677 static PetscErrorCode MatCreateSubMatrices_SeqDense(Mat A,PetscInt n,const IS irow[],const IS icol[],MatReuse scall,Mat *B[]) 1678 { 1679 PetscErrorCode ierr; 1680 PetscInt i; 1681 1682 PetscFunctionBegin; 1683 if (scall == MAT_INITIAL_MATRIX) { 1684 ierr = PetscCalloc1(n+1,B);CHKERRQ(ierr); 1685 } 1686 1687 for (i=0; i<n; i++) { 1688 ierr = MatCreateSubMatrix_SeqDense(A,irow[i],icol[i],PETSC_DECIDE,scall,&(*B)[i]);CHKERRQ(ierr); 1689 } 1690 PetscFunctionReturn(0); 1691 } 1692 1693 static PetscErrorCode MatAssemblyBegin_SeqDense(Mat mat,MatAssemblyType mode) 1694 { 1695 PetscFunctionBegin; 1696 PetscFunctionReturn(0); 1697 } 1698 1699 static PetscErrorCode MatAssemblyEnd_SeqDense(Mat mat,MatAssemblyType mode) 1700 { 1701 PetscFunctionBegin; 1702 PetscFunctionReturn(0); 1703 } 1704 1705 static PetscErrorCode MatCopy_SeqDense(Mat A,Mat B,MatStructure str) 1706 { 1707 Mat_SeqDense *a = (Mat_SeqDense*)A->data,*b = (Mat_SeqDense*)B->data; 1708 PetscErrorCode ierr; 1709 PetscInt lda1=a->lda,lda2=b->lda, m=A->rmap->n,n=A->cmap->n, j; 1710 1711 PetscFunctionBegin; 1712 /* If the two matrices don't have the same copy implementation, they aren't compatible for fast copy. */ 1713 if (A->ops->copy != B->ops->copy) { 1714 ierr = MatCopy_Basic(A,B,str);CHKERRQ(ierr); 1715 PetscFunctionReturn(0); 1716 } 1717 if (m != B->rmap->n || n != B->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"size(B) != size(A)"); 1718 if (lda1>m || lda2>m) { 1719 for (j=0; j<n; j++) { 1720 ierr = PetscMemcpy(b->v+j*lda2,a->v+j*lda1,m*sizeof(PetscScalar));CHKERRQ(ierr); 1721 } 1722 } else { 1723 ierr = PetscMemcpy(b->v,a->v,A->rmap->n*A->cmap->n*sizeof(PetscScalar));CHKERRQ(ierr); 1724 } 1725 PetscFunctionReturn(0); 1726 } 1727 1728 static PetscErrorCode MatSetUp_SeqDense(Mat A) 1729 { 1730 PetscErrorCode ierr; 1731 1732 PetscFunctionBegin; 1733 ierr = MatSeqDenseSetPreallocation(A,0);CHKERRQ(ierr); 1734 PetscFunctionReturn(0); 1735 } 1736 1737 static PetscErrorCode MatConjugate_SeqDense(Mat A) 1738 { 1739 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1740 PetscInt i,nz = A->rmap->n*A->cmap->n; 1741 PetscScalar *aa = a->v; 1742 1743 PetscFunctionBegin; 1744 for (i=0; i<nz; i++) aa[i] = PetscConj(aa[i]); 1745 PetscFunctionReturn(0); 1746 } 1747 1748 static PetscErrorCode MatRealPart_SeqDense(Mat A) 1749 { 1750 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1751 PetscInt i,nz = A->rmap->n*A->cmap->n; 1752 PetscScalar *aa = a->v; 1753 1754 PetscFunctionBegin; 1755 for (i=0; i<nz; i++) aa[i] = PetscRealPart(aa[i]); 1756 PetscFunctionReturn(0); 1757 } 1758 1759 static PetscErrorCode MatImaginaryPart_SeqDense(Mat A) 1760 { 1761 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1762 PetscInt i,nz = A->rmap->n*A->cmap->n; 1763 PetscScalar *aa = a->v; 1764 1765 PetscFunctionBegin; 1766 for (i=0; i<nz; i++) aa[i] = PetscImaginaryPart(aa[i]); 1767 PetscFunctionReturn(0); 1768 } 1769 1770 /* ----------------------------------------------------------------*/ 1771 PETSC_INTERN PetscErrorCode MatMatMult_SeqDense_SeqDense(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C) 1772 { 1773 PetscErrorCode ierr; 1774 1775 PetscFunctionBegin; 1776 if (scall == MAT_INITIAL_MATRIX) { 1777 ierr = PetscLogEventBegin(MAT_MatMultSymbolic,A,B,0,0);CHKERRQ(ierr); 1778 ierr = MatMatMultSymbolic_SeqDense_SeqDense(A,B,fill,C);CHKERRQ(ierr); 1779 ierr = PetscLogEventEnd(MAT_MatMultSymbolic,A,B,0,0);CHKERRQ(ierr); 1780 } 1781 ierr = PetscLogEventBegin(MAT_MatMultNumeric,A,B,0,0);CHKERRQ(ierr); 1782 ierr = MatMatMultNumeric_SeqDense_SeqDense(A,B,*C);CHKERRQ(ierr); 1783 ierr = PetscLogEventEnd(MAT_MatMultNumeric,A,B,0,0);CHKERRQ(ierr); 1784 PetscFunctionReturn(0); 1785 } 1786 1787 PetscErrorCode MatMatMultSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat *C) 1788 { 1789 PetscErrorCode ierr; 1790 PetscInt m=A->rmap->n,n=B->cmap->n; 1791 Mat Cmat; 1792 1793 PetscFunctionBegin; 1794 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); 1795 ierr = MatCreate(PETSC_COMM_SELF,&Cmat);CHKERRQ(ierr); 1796 ierr = MatSetSizes(Cmat,m,n,m,n);CHKERRQ(ierr); 1797 ierr = MatSetType(Cmat,MATSEQDENSE);CHKERRQ(ierr); 1798 ierr = MatSeqDenseSetPreallocation(Cmat,NULL);CHKERRQ(ierr); 1799 1800 *C = Cmat; 1801 PetscFunctionReturn(0); 1802 } 1803 1804 PetscErrorCode MatMatMultNumeric_SeqDense_SeqDense(Mat A,Mat B,Mat C) 1805 { 1806 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1807 Mat_SeqDense *b = (Mat_SeqDense*)B->data; 1808 Mat_SeqDense *c = (Mat_SeqDense*)C->data; 1809 PetscBLASInt m,n,k; 1810 PetscScalar _DOne=1.0,_DZero=0.0; 1811 PetscErrorCode ierr; 1812 PetscBool flg; 1813 1814 PetscFunctionBegin; 1815 ierr = PetscObjectTypeCompare((PetscObject)B,MATSEQDENSE,&flg);CHKERRQ(ierr); 1816 if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Second matrix must be dense"); 1817 1818 /* Handle case where where user provided the final C matrix rather than calling MatMatMult() with MAT_INITIAL_MATRIX*/ 1819 ierr = PetscObjectTypeCompare((PetscObject)A,MATSEQAIJ,&flg);CHKERRQ(ierr); 1820 if (flg) { 1821 C->ops->matmultnumeric = MatMatMultNumeric_SeqAIJ_SeqDense; 1822 ierr = (*C->ops->matmultnumeric)(A,B,C);CHKERRQ(ierr); 1823 PetscFunctionReturn(0); 1824 } 1825 1826 ierr = PetscObjectTypeCompare((PetscObject)A,MATSEQDENSE,&flg);CHKERRQ(ierr); 1827 if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"First matrix must be dense"); 1828 ierr = PetscBLASIntCast(A->rmap->n,&m);CHKERRQ(ierr); 1829 ierr = PetscBLASIntCast(B->cmap->n,&n);CHKERRQ(ierr); 1830 ierr = PetscBLASIntCast(A->cmap->n,&k);CHKERRQ(ierr); 1831 PetscStackCallBLAS("BLASgemm",BLASgemm_("N","N",&m,&n,&k,&_DOne,a->v,&a->lda,b->v,&b->lda,&_DZero,c->v,&c->lda)); 1832 PetscFunctionReturn(0); 1833 } 1834 1835 PetscErrorCode MatTransposeMatMult_SeqDense_SeqDense(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C) 1836 { 1837 PetscErrorCode ierr; 1838 1839 PetscFunctionBegin; 1840 if (scall == MAT_INITIAL_MATRIX) { 1841 ierr = PetscLogEventBegin(MAT_TransposeMatMultSymbolic,A,B,0,0);CHKERRQ(ierr); 1842 ierr = MatTransposeMatMultSymbolic_SeqDense_SeqDense(A,B,fill,C);CHKERRQ(ierr); 1843 ierr = PetscLogEventEnd(MAT_TransposeMatMultSymbolic,A,B,0,0);CHKERRQ(ierr); 1844 } 1845 ierr = PetscLogEventBegin(MAT_TransposeMatMultNumeric,A,B,0,0);CHKERRQ(ierr); 1846 ierr = MatTransposeMatMultNumeric_SeqDense_SeqDense(A,B,*C);CHKERRQ(ierr); 1847 ierr = PetscLogEventEnd(MAT_TransposeMatMultNumeric,A,B,0,0);CHKERRQ(ierr); 1848 PetscFunctionReturn(0); 1849 } 1850 1851 PetscErrorCode MatTransposeMatMultSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat *C) 1852 { 1853 PetscErrorCode ierr; 1854 PetscInt m=A->cmap->n,n=B->cmap->n; 1855 Mat Cmat; 1856 1857 PetscFunctionBegin; 1858 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); 1859 ierr = MatCreate(PETSC_COMM_SELF,&Cmat);CHKERRQ(ierr); 1860 ierr = MatSetSizes(Cmat,m,n,m,n);CHKERRQ(ierr); 1861 ierr = MatSetType(Cmat,MATSEQDENSE);CHKERRQ(ierr); 1862 ierr = MatSeqDenseSetPreallocation(Cmat,NULL);CHKERRQ(ierr); 1863 1864 Cmat->assembled = PETSC_TRUE; 1865 1866 *C = Cmat; 1867 PetscFunctionReturn(0); 1868 } 1869 1870 PetscErrorCode MatTransposeMatMultNumeric_SeqDense_SeqDense(Mat A,Mat B,Mat C) 1871 { 1872 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1873 Mat_SeqDense *b = (Mat_SeqDense*)B->data; 1874 Mat_SeqDense *c = (Mat_SeqDense*)C->data; 1875 PetscBLASInt m,n,k; 1876 PetscScalar _DOne=1.0,_DZero=0.0; 1877 PetscErrorCode ierr; 1878 1879 PetscFunctionBegin; 1880 ierr = PetscBLASIntCast(A->cmap->n,&m);CHKERRQ(ierr); 1881 ierr = PetscBLASIntCast(B->cmap->n,&n);CHKERRQ(ierr); 1882 ierr = PetscBLASIntCast(A->rmap->n,&k);CHKERRQ(ierr); 1883 /* 1884 Note the m and n arguments below are the number rows and columns of A', not A! 1885 */ 1886 PetscStackCallBLAS("BLASgemm",BLASgemm_("T","N",&m,&n,&k,&_DOne,a->v,&a->lda,b->v,&b->lda,&_DZero,c->v,&c->lda)); 1887 PetscFunctionReturn(0); 1888 } 1889 1890 static PetscErrorCode MatGetRowMax_SeqDense(Mat A,Vec v,PetscInt idx[]) 1891 { 1892 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1893 PetscErrorCode ierr; 1894 PetscInt i,j,m = A->rmap->n,n = A->cmap->n,p; 1895 PetscScalar *x; 1896 MatScalar *aa = a->v; 1897 1898 PetscFunctionBegin; 1899 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 1900 1901 ierr = VecSet(v,0.0);CHKERRQ(ierr); 1902 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 1903 ierr = VecGetLocalSize(v,&p);CHKERRQ(ierr); 1904 if (p != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 1905 for (i=0; i<m; i++) { 1906 x[i] = aa[i]; if (idx) idx[i] = 0; 1907 for (j=1; j<n; j++) { 1908 if (PetscRealPart(x[i]) < PetscRealPart(aa[i+m*j])) {x[i] = aa[i + m*j]; if (idx) idx[i] = j;} 1909 } 1910 } 1911 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 1912 PetscFunctionReturn(0); 1913 } 1914 1915 static PetscErrorCode MatGetRowMaxAbs_SeqDense(Mat A,Vec v,PetscInt idx[]) 1916 { 1917 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1918 PetscErrorCode ierr; 1919 PetscInt i,j,m = A->rmap->n,n = A->cmap->n,p; 1920 PetscScalar *x; 1921 PetscReal atmp; 1922 MatScalar *aa = a->v; 1923 1924 PetscFunctionBegin; 1925 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 1926 1927 ierr = VecSet(v,0.0);CHKERRQ(ierr); 1928 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 1929 ierr = VecGetLocalSize(v,&p);CHKERRQ(ierr); 1930 if (p != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 1931 for (i=0; i<m; i++) { 1932 x[i] = PetscAbsScalar(aa[i]); 1933 for (j=1; j<n; j++) { 1934 atmp = PetscAbsScalar(aa[i+m*j]); 1935 if (PetscAbsScalar(x[i]) < atmp) {x[i] = atmp; if (idx) idx[i] = j;} 1936 } 1937 } 1938 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 1939 PetscFunctionReturn(0); 1940 } 1941 1942 static PetscErrorCode MatGetRowMin_SeqDense(Mat A,Vec v,PetscInt idx[]) 1943 { 1944 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1945 PetscErrorCode ierr; 1946 PetscInt i,j,m = A->rmap->n,n = A->cmap->n,p; 1947 PetscScalar *x; 1948 MatScalar *aa = a->v; 1949 1950 PetscFunctionBegin; 1951 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 1952 1953 ierr = VecSet(v,0.0);CHKERRQ(ierr); 1954 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 1955 ierr = VecGetLocalSize(v,&p);CHKERRQ(ierr); 1956 if (p != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 1957 for (i=0; i<m; i++) { 1958 x[i] = aa[i]; if (idx) idx[i] = 0; 1959 for (j=1; j<n; j++) { 1960 if (PetscRealPart(x[i]) > PetscRealPart(aa[i+m*j])) {x[i] = aa[i + m*j]; if (idx) idx[i] = j;} 1961 } 1962 } 1963 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 1964 PetscFunctionReturn(0); 1965 } 1966 1967 static PetscErrorCode MatGetColumnVector_SeqDense(Mat A,Vec v,PetscInt col) 1968 { 1969 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1970 PetscErrorCode ierr; 1971 PetscScalar *x; 1972 1973 PetscFunctionBegin; 1974 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 1975 1976 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 1977 ierr = PetscMemcpy(x,a->v+col*a->lda,A->rmap->n*sizeof(PetscScalar));CHKERRQ(ierr); 1978 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 1979 PetscFunctionReturn(0); 1980 } 1981 1982 1983 PetscErrorCode MatGetColumnNorms_SeqDense(Mat A,NormType type,PetscReal *norms) 1984 { 1985 PetscErrorCode ierr; 1986 PetscInt i,j,m,n; 1987 PetscScalar *a; 1988 1989 PetscFunctionBegin; 1990 ierr = MatGetSize(A,&m,&n);CHKERRQ(ierr); 1991 ierr = PetscMemzero(norms,n*sizeof(PetscReal));CHKERRQ(ierr); 1992 ierr = MatDenseGetArray(A,&a);CHKERRQ(ierr); 1993 if (type == NORM_2) { 1994 for (i=0; i<n; i++) { 1995 for (j=0; j<m; j++) { 1996 norms[i] += PetscAbsScalar(a[j]*a[j]); 1997 } 1998 a += m; 1999 } 2000 } else if (type == NORM_1) { 2001 for (i=0; i<n; i++) { 2002 for (j=0; j<m; j++) { 2003 norms[i] += PetscAbsScalar(a[j]); 2004 } 2005 a += m; 2006 } 2007 } else if (type == NORM_INFINITY) { 2008 for (i=0; i<n; i++) { 2009 for (j=0; j<m; j++) { 2010 norms[i] = PetscMax(PetscAbsScalar(a[j]),norms[i]); 2011 } 2012 a += m; 2013 } 2014 } else SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Unknown NormType"); 2015 ierr = MatDenseRestoreArray(A,&a);CHKERRQ(ierr); 2016 if (type == NORM_2) { 2017 for (i=0; i<n; i++) norms[i] = PetscSqrtReal(norms[i]); 2018 } 2019 PetscFunctionReturn(0); 2020 } 2021 2022 static PetscErrorCode MatSetRandom_SeqDense(Mat x,PetscRandom rctx) 2023 { 2024 PetscErrorCode ierr; 2025 PetscScalar *a; 2026 PetscInt m,n,i; 2027 2028 PetscFunctionBegin; 2029 ierr = MatGetSize(x,&m,&n);CHKERRQ(ierr); 2030 ierr = MatDenseGetArray(x,&a);CHKERRQ(ierr); 2031 for (i=0; i<m*n; i++) { 2032 ierr = PetscRandomGetValue(rctx,a+i);CHKERRQ(ierr); 2033 } 2034 ierr = MatDenseRestoreArray(x,&a);CHKERRQ(ierr); 2035 PetscFunctionReturn(0); 2036 } 2037 2038 static PetscErrorCode MatMissingDiagonal_SeqDense(Mat A,PetscBool *missing,PetscInt *d) 2039 { 2040 PetscFunctionBegin; 2041 *missing = PETSC_FALSE; 2042 PetscFunctionReturn(0); 2043 } 2044 2045 2046 /* -------------------------------------------------------------------*/ 2047 static struct _MatOps MatOps_Values = { MatSetValues_SeqDense, 2048 MatGetRow_SeqDense, 2049 MatRestoreRow_SeqDense, 2050 MatMult_SeqDense, 2051 /* 4*/ MatMultAdd_SeqDense, 2052 MatMultTranspose_SeqDense, 2053 MatMultTransposeAdd_SeqDense, 2054 0, 2055 0, 2056 0, 2057 /* 10*/ 0, 2058 MatLUFactor_SeqDense, 2059 MatCholeskyFactor_SeqDense, 2060 MatSOR_SeqDense, 2061 MatTranspose_SeqDense, 2062 /* 15*/ MatGetInfo_SeqDense, 2063 MatEqual_SeqDense, 2064 MatGetDiagonal_SeqDense, 2065 MatDiagonalScale_SeqDense, 2066 MatNorm_SeqDense, 2067 /* 20*/ MatAssemblyBegin_SeqDense, 2068 MatAssemblyEnd_SeqDense, 2069 MatSetOption_SeqDense, 2070 MatZeroEntries_SeqDense, 2071 /* 24*/ MatZeroRows_SeqDense, 2072 0, 2073 0, 2074 0, 2075 0, 2076 /* 29*/ MatSetUp_SeqDense, 2077 0, 2078 0, 2079 0, 2080 0, 2081 /* 34*/ MatDuplicate_SeqDense, 2082 0, 2083 0, 2084 0, 2085 0, 2086 /* 39*/ MatAXPY_SeqDense, 2087 MatCreateSubMatrices_SeqDense, 2088 0, 2089 MatGetValues_SeqDense, 2090 MatCopy_SeqDense, 2091 /* 44*/ MatGetRowMax_SeqDense, 2092 MatScale_SeqDense, 2093 MatShift_Basic, 2094 0, 2095 MatZeroRowsColumns_SeqDense, 2096 /* 49*/ MatSetRandom_SeqDense, 2097 0, 2098 0, 2099 0, 2100 0, 2101 /* 54*/ 0, 2102 0, 2103 0, 2104 0, 2105 0, 2106 /* 59*/ 0, 2107 MatDestroy_SeqDense, 2108 MatView_SeqDense, 2109 0, 2110 0, 2111 /* 64*/ 0, 2112 0, 2113 0, 2114 0, 2115 0, 2116 /* 69*/ MatGetRowMaxAbs_SeqDense, 2117 0, 2118 0, 2119 0, 2120 0, 2121 /* 74*/ 0, 2122 0, 2123 0, 2124 0, 2125 0, 2126 /* 79*/ 0, 2127 0, 2128 0, 2129 0, 2130 /* 83*/ MatLoad_SeqDense, 2131 0, 2132 MatIsHermitian_SeqDense, 2133 0, 2134 0, 2135 0, 2136 /* 89*/ MatMatMult_SeqDense_SeqDense, 2137 MatMatMultSymbolic_SeqDense_SeqDense, 2138 MatMatMultNumeric_SeqDense_SeqDense, 2139 MatPtAP_SeqDense_SeqDense, 2140 MatPtAPSymbolic_SeqDense_SeqDense, 2141 /* 94*/ MatPtAPNumeric_SeqDense_SeqDense, 2142 0, 2143 0, 2144 0, 2145 0, 2146 /* 99*/ 0, 2147 0, 2148 0, 2149 MatConjugate_SeqDense, 2150 0, 2151 /*104*/ 0, 2152 MatRealPart_SeqDense, 2153 MatImaginaryPart_SeqDense, 2154 0, 2155 0, 2156 /*109*/ MatMatSolve_SeqDense, 2157 0, 2158 MatGetRowMin_SeqDense, 2159 MatGetColumnVector_SeqDense, 2160 MatMissingDiagonal_SeqDense, 2161 /*114*/ 0, 2162 0, 2163 0, 2164 0, 2165 0, 2166 /*119*/ 0, 2167 0, 2168 0, 2169 0, 2170 0, 2171 /*124*/ 0, 2172 MatGetColumnNorms_SeqDense, 2173 0, 2174 0, 2175 0, 2176 /*129*/ 0, 2177 MatTransposeMatMult_SeqDense_SeqDense, 2178 MatTransposeMatMultSymbolic_SeqDense_SeqDense, 2179 MatTransposeMatMultNumeric_SeqDense_SeqDense, 2180 0, 2181 /*134*/ 0, 2182 0, 2183 0, 2184 0, 2185 0, 2186 /*139*/ 0, 2187 0, 2188 0 2189 }; 2190 2191 /*@C 2192 MatCreateSeqDense - Creates a sequential dense matrix that 2193 is stored in column major order (the usual Fortran 77 manner). Many 2194 of the matrix operations use the BLAS and LAPACK routines. 2195 2196 Collective on MPI_Comm 2197 2198 Input Parameters: 2199 + comm - MPI communicator, set to PETSC_COMM_SELF 2200 . m - number of rows 2201 . n - number of columns 2202 - data - optional location of matrix data in column major order. Set data=NULL for PETSc 2203 to control all matrix memory allocation. 2204 2205 Output Parameter: 2206 . A - the matrix 2207 2208 Notes: 2209 The data input variable is intended primarily for Fortran programmers 2210 who wish to allocate their own matrix memory space. Most users should 2211 set data=NULL. 2212 2213 Level: intermediate 2214 2215 .keywords: dense, matrix, LAPACK, BLAS 2216 2217 .seealso: MatCreate(), MatCreateDense(), MatSetValues() 2218 @*/ 2219 PetscErrorCode MatCreateSeqDense(MPI_Comm comm,PetscInt m,PetscInt n,PetscScalar *data,Mat *A) 2220 { 2221 PetscErrorCode ierr; 2222 2223 PetscFunctionBegin; 2224 ierr = MatCreate(comm,A);CHKERRQ(ierr); 2225 ierr = MatSetSizes(*A,m,n,m,n);CHKERRQ(ierr); 2226 ierr = MatSetType(*A,MATSEQDENSE);CHKERRQ(ierr); 2227 ierr = MatSeqDenseSetPreallocation(*A,data);CHKERRQ(ierr); 2228 PetscFunctionReturn(0); 2229 } 2230 2231 /*@C 2232 MatSeqDenseSetPreallocation - Sets the array used for storing the matrix elements 2233 2234 Collective on MPI_Comm 2235 2236 Input Parameters: 2237 + B - the matrix 2238 - data - the array (or NULL) 2239 2240 Notes: 2241 The data input variable is intended primarily for Fortran programmers 2242 who wish to allocate their own matrix memory space. Most users should 2243 need not call this routine. 2244 2245 Level: intermediate 2246 2247 .keywords: dense, matrix, LAPACK, BLAS 2248 2249 .seealso: MatCreate(), MatCreateDense(), MatSetValues(), MatSeqDenseSetLDA() 2250 2251 @*/ 2252 PetscErrorCode MatSeqDenseSetPreallocation(Mat B,PetscScalar data[]) 2253 { 2254 PetscErrorCode ierr; 2255 2256 PetscFunctionBegin; 2257 ierr = PetscTryMethod(B,"MatSeqDenseSetPreallocation_C",(Mat,PetscScalar[]),(B,data));CHKERRQ(ierr); 2258 PetscFunctionReturn(0); 2259 } 2260 2261 PetscErrorCode MatSeqDenseSetPreallocation_SeqDense(Mat B,PetscScalar *data) 2262 { 2263 Mat_SeqDense *b; 2264 PetscErrorCode ierr; 2265 2266 PetscFunctionBegin; 2267 B->preallocated = PETSC_TRUE; 2268 2269 ierr = PetscLayoutSetUp(B->rmap);CHKERRQ(ierr); 2270 ierr = PetscLayoutSetUp(B->cmap);CHKERRQ(ierr); 2271 2272 b = (Mat_SeqDense*)B->data; 2273 b->Mmax = B->rmap->n; 2274 b->Nmax = B->cmap->n; 2275 if (b->lda <= 0 || b->changelda) b->lda = B->rmap->n; 2276 2277 ierr = PetscIntMultError(b->lda,b->Nmax,NULL);CHKERRQ(ierr); 2278 if (!data) { /* petsc-allocated storage */ 2279 if (!b->user_alloc) { ierr = PetscFree(b->v);CHKERRQ(ierr); } 2280 ierr = PetscCalloc1((size_t)b->lda*b->Nmax,&b->v);CHKERRQ(ierr); 2281 ierr = PetscLogObjectMemory((PetscObject)B,b->lda*b->Nmax*sizeof(PetscScalar));CHKERRQ(ierr); 2282 2283 b->user_alloc = PETSC_FALSE; 2284 } else { /* user-allocated storage */ 2285 if (!b->user_alloc) { ierr = PetscFree(b->v);CHKERRQ(ierr); } 2286 b->v = data; 2287 b->user_alloc = PETSC_TRUE; 2288 } 2289 B->assembled = PETSC_TRUE; 2290 PetscFunctionReturn(0); 2291 } 2292 2293 #if defined(PETSC_HAVE_ELEMENTAL) 2294 PETSC_INTERN PetscErrorCode MatConvert_SeqDense_Elemental(Mat A, MatType newtype,MatReuse reuse,Mat *newmat) 2295 { 2296 Mat mat_elemental; 2297 PetscErrorCode ierr; 2298 PetscScalar *array,*v_colwise; 2299 PetscInt M=A->rmap->N,N=A->cmap->N,i,j,k,*rows,*cols; 2300 2301 PetscFunctionBegin; 2302 ierr = PetscMalloc3(M*N,&v_colwise,M,&rows,N,&cols);CHKERRQ(ierr); 2303 ierr = MatDenseGetArray(A,&array);CHKERRQ(ierr); 2304 /* convert column-wise array into row-wise v_colwise, see MatSetValues_Elemental() */ 2305 k = 0; 2306 for (j=0; j<N; j++) { 2307 cols[j] = j; 2308 for (i=0; i<M; i++) { 2309 v_colwise[j*M+i] = array[k++]; 2310 } 2311 } 2312 for (i=0; i<M; i++) { 2313 rows[i] = i; 2314 } 2315 ierr = MatDenseRestoreArray(A,&array);CHKERRQ(ierr); 2316 2317 ierr = MatCreate(PetscObjectComm((PetscObject)A), &mat_elemental);CHKERRQ(ierr); 2318 ierr = MatSetSizes(mat_elemental,PETSC_DECIDE,PETSC_DECIDE,M,N);CHKERRQ(ierr); 2319 ierr = MatSetType(mat_elemental,MATELEMENTAL);CHKERRQ(ierr); 2320 ierr = MatSetUp(mat_elemental);CHKERRQ(ierr); 2321 2322 /* PETSc-Elemental interaface uses axpy for setting off-processor entries, only ADD_VALUES is allowed */ 2323 ierr = MatSetValues(mat_elemental,M,rows,N,cols,v_colwise,ADD_VALUES);CHKERRQ(ierr); 2324 ierr = MatAssemblyBegin(mat_elemental, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2325 ierr = MatAssemblyEnd(mat_elemental, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2326 ierr = PetscFree3(v_colwise,rows,cols);CHKERRQ(ierr); 2327 2328 if (reuse == MAT_INPLACE_MATRIX) { 2329 ierr = MatHeaderReplace(A,&mat_elemental);CHKERRQ(ierr); 2330 } else { 2331 *newmat = mat_elemental; 2332 } 2333 PetscFunctionReturn(0); 2334 } 2335 #endif 2336 2337 /*@C 2338 MatSeqDenseSetLDA - Declare the leading dimension of the user-provided array 2339 2340 Input parameter: 2341 + A - the matrix 2342 - lda - the leading dimension 2343 2344 Notes: 2345 This routine is to be used in conjunction with MatSeqDenseSetPreallocation(); 2346 it asserts that the preallocation has a leading dimension (the LDA parameter 2347 of Blas and Lapack fame) larger than M, the first dimension of the matrix. 2348 2349 Level: intermediate 2350 2351 .keywords: dense, matrix, LAPACK, BLAS 2352 2353 .seealso: MatCreate(), MatCreateSeqDense(), MatSeqDenseSetPreallocation(), MatSetMaximumSize() 2354 2355 @*/ 2356 PetscErrorCode MatSeqDenseSetLDA(Mat B,PetscInt lda) 2357 { 2358 Mat_SeqDense *b = (Mat_SeqDense*)B->data; 2359 2360 PetscFunctionBegin; 2361 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); 2362 b->lda = lda; 2363 b->changelda = PETSC_FALSE; 2364 b->Mmax = PetscMax(b->Mmax,lda); 2365 PetscFunctionReturn(0); 2366 } 2367 2368 /*MC 2369 MATSEQDENSE - MATSEQDENSE = "seqdense" - A matrix type to be used for sequential dense matrices. 2370 2371 Options Database Keys: 2372 . -mat_type seqdense - sets the matrix type to "seqdense" during a call to MatSetFromOptions() 2373 2374 Level: beginner 2375 2376 .seealso: MatCreateSeqDense() 2377 2378 M*/ 2379 2380 PETSC_EXTERN PetscErrorCode MatCreate_SeqDense(Mat B) 2381 { 2382 Mat_SeqDense *b; 2383 PetscErrorCode ierr; 2384 PetscMPIInt size; 2385 2386 PetscFunctionBegin; 2387 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)B),&size);CHKERRQ(ierr); 2388 if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Comm must be of size 1"); 2389 2390 ierr = PetscNewLog(B,&b);CHKERRQ(ierr); 2391 ierr = PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 2392 B->data = (void*)b; 2393 2394 b->pivots = 0; 2395 b->roworiented = PETSC_TRUE; 2396 b->v = 0; 2397 b->changelda = PETSC_FALSE; 2398 2399 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseGetArray_C",MatDenseGetArray_SeqDense);CHKERRQ(ierr); 2400 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDensePlaceArray_C",MatDensePlaceArray_SeqDense);CHKERRQ(ierr); 2401 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseResetArray_C",MatDenseResetArray_SeqDense);CHKERRQ(ierr); 2402 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreArray_C",MatDenseRestoreArray_SeqDense);CHKERRQ(ierr); 2403 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqdense_seqaij_C",MatConvert_SeqDense_SeqAIJ);CHKERRQ(ierr); 2404 #if defined(PETSC_HAVE_ELEMENTAL) 2405 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqdense_elemental_C",MatConvert_SeqDense_Elemental);CHKERRQ(ierr); 2406 #endif 2407 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqDenseSetPreallocation_C",MatSeqDenseSetPreallocation_SeqDense);CHKERRQ(ierr); 2408 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMult_seqaij_seqdense_C",MatMatMult_SeqAIJ_SeqDense);CHKERRQ(ierr); 2409 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMultSymbolic_seqaij_seqdense_C",MatMatMultSymbolic_SeqAIJ_SeqDense);CHKERRQ(ierr); 2410 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMultNumeric_seqaij_seqdense_C",MatMatMultNumeric_SeqAIJ_SeqDense);CHKERRQ(ierr); 2411 ierr = PetscObjectComposeFunction((PetscObject)B,"MatPtAP_seqaij_seqdense_C",MatPtAP_SeqDense_SeqDense);CHKERRQ(ierr); 2412 ierr = PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMult_seqaij_seqdense_C",MatTransposeMatMult_SeqAIJ_SeqDense);CHKERRQ(ierr); 2413 ierr = PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultSymbolic_seqaij_seqdense_C",MatTransposeMatMultSymbolic_SeqAIJ_SeqDense);CHKERRQ(ierr); 2414 ierr = PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultNumeric_seqaij_seqdense_C",MatTransposeMatMultNumeric_SeqAIJ_SeqDense);CHKERRQ(ierr); 2415 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQDENSE);CHKERRQ(ierr); 2416 PetscFunctionReturn(0); 2417 } 2418