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