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