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