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