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