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