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