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