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