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