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 ierr = PetscBLASIntCast(m,&bm);CHKERRQ(ierr); 600 if (flag & SOR_ZERO_INITIAL_GUESS) { 601 /* this is a hack fix, should have another version without the second BLASdot */ 602 ierr = VecSet(xx,zero);CHKERRQ(ierr); 603 } 604 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 605 ierr = VecGetArrayRead(bb,&b);CHKERRQ(ierr); 606 its = its*lits; 607 if (its <= 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Relaxation requires global its %D and local its %D both positive",its,lits); 608 while (its--) { 609 if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP) { 610 for (i=0; i<m; i++) { 611 PetscStackCallBLAS("BLASdot",xt = b[i] - BLASdot_(&bm,v+i,&bm,x,&o)); 612 x[i] = (1. - omega)*x[i] + omega*(xt+v[i + i*m]*x[i])/(v[i + i*m]+shift); 613 } 614 } 615 if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP) { 616 for (i=m-1; i>=0; i--) { 617 PetscStackCallBLAS("BLASdot",xt = b[i] - BLASdot_(&bm,v+i,&bm,x,&o)); 618 x[i] = (1. - omega)*x[i] + omega*(xt+v[i + i*m]*x[i])/(v[i + i*m]+shift); 619 } 620 } 621 } 622 ierr = VecRestoreArrayRead(bb,&b);CHKERRQ(ierr); 623 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 624 PetscFunctionReturn(0); 625 } 626 627 /* -----------------------------------------------------------------*/ 628 #undef __FUNCT__ 629 #define __FUNCT__ "MatMultTranspose_SeqDense" 630 PetscErrorCode MatMultTranspose_SeqDense(Mat A,Vec xx,Vec yy) 631 { 632 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 633 const PetscScalar *v = mat->v,*x; 634 PetscScalar *y; 635 PetscErrorCode ierr; 636 PetscBLASInt m, n,_One=1; 637 PetscScalar _DOne=1.0,_DZero=0.0; 638 639 PetscFunctionBegin; 640 ierr = PetscBLASIntCast(A->rmap->n,&m);CHKERRQ(ierr); 641 ierr = PetscBLASIntCast(A->cmap->n,&n);CHKERRQ(ierr); 642 if (!A->rmap->n || !A->cmap->n) PetscFunctionReturn(0); 643 ierr = VecGetArrayRead(xx,&x);CHKERRQ(ierr); 644 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 645 PetscStackCallBLAS("BLASgemv",BLASgemv_("T",&m,&n,&_DOne,v,&mat->lda,x,&_One,&_DZero,y,&_One)); 646 ierr = VecRestoreArrayRead(xx,&x);CHKERRQ(ierr); 647 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 648 ierr = PetscLogFlops(2.0*A->rmap->n*A->cmap->n - A->cmap->n);CHKERRQ(ierr); 649 PetscFunctionReturn(0); 650 } 651 652 #undef __FUNCT__ 653 #define __FUNCT__ "MatMult_SeqDense" 654 PetscErrorCode MatMult_SeqDense(Mat A,Vec xx,Vec yy) 655 { 656 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 657 PetscScalar *y,_DOne=1.0,_DZero=0.0; 658 PetscErrorCode ierr; 659 PetscBLASInt m, n, _One=1; 660 const PetscScalar *v = mat->v,*x; 661 662 PetscFunctionBegin; 663 ierr = PetscBLASIntCast(A->rmap->n,&m);CHKERRQ(ierr); 664 ierr = PetscBLASIntCast(A->cmap->n,&n);CHKERRQ(ierr); 665 if (!A->rmap->n || !A->cmap->n) PetscFunctionReturn(0); 666 ierr = VecGetArrayRead(xx,&x);CHKERRQ(ierr); 667 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 668 PetscStackCallBLAS("BLASgemv",BLASgemv_("N",&m,&n,&_DOne,v,&(mat->lda),x,&_One,&_DZero,y,&_One)); 669 ierr = VecRestoreArrayRead(xx,&x);CHKERRQ(ierr); 670 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 671 ierr = PetscLogFlops(2.0*A->rmap->n*A->cmap->n - A->rmap->n);CHKERRQ(ierr); 672 PetscFunctionReturn(0); 673 } 674 675 #undef __FUNCT__ 676 #define __FUNCT__ "MatMultAdd_SeqDense" 677 PetscErrorCode MatMultAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy) 678 { 679 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 680 const PetscScalar *v = mat->v,*x; 681 PetscScalar *y,_DOne=1.0; 682 PetscErrorCode ierr; 683 PetscBLASInt m, n, _One=1; 684 685 PetscFunctionBegin; 686 ierr = PetscBLASIntCast(A->rmap->n,&m);CHKERRQ(ierr); 687 ierr = PetscBLASIntCast(A->cmap->n,&n);CHKERRQ(ierr); 688 if (!A->rmap->n || !A->cmap->n) PetscFunctionReturn(0); 689 if (zz != yy) {ierr = VecCopy(zz,yy);CHKERRQ(ierr);} 690 ierr = VecGetArrayRead(xx,&x);CHKERRQ(ierr); 691 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 692 PetscStackCallBLAS("BLASgemv",BLASgemv_("N",&m,&n,&_DOne,v,&(mat->lda),x,&_One,&_DOne,y,&_One)); 693 ierr = VecRestoreArrayRead(xx,&x);CHKERRQ(ierr); 694 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 695 ierr = PetscLogFlops(2.0*A->rmap->n*A->cmap->n);CHKERRQ(ierr); 696 PetscFunctionReturn(0); 697 } 698 699 #undef __FUNCT__ 700 #define __FUNCT__ "MatMultTransposeAdd_SeqDense" 701 PetscErrorCode MatMultTransposeAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy) 702 { 703 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 704 const PetscScalar *v = mat->v,*x; 705 PetscScalar *y; 706 PetscErrorCode ierr; 707 PetscBLASInt m, n, _One=1; 708 PetscScalar _DOne=1.0; 709 710 PetscFunctionBegin; 711 ierr = PetscBLASIntCast(A->rmap->n,&m);CHKERRQ(ierr); 712 ierr = PetscBLASIntCast(A->cmap->n,&n);CHKERRQ(ierr); 713 if (!A->rmap->n || !A->cmap->n) PetscFunctionReturn(0); 714 if (zz != yy) {ierr = VecCopy(zz,yy);CHKERRQ(ierr);} 715 ierr = VecGetArrayRead(xx,&x);CHKERRQ(ierr); 716 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 717 PetscStackCallBLAS("BLASgemv",BLASgemv_("T",&m,&n,&_DOne,v,&(mat->lda),x,&_One,&_DOne,y,&_One)); 718 ierr = VecRestoreArrayRead(xx,&x);CHKERRQ(ierr); 719 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 720 ierr = PetscLogFlops(2.0*A->rmap->n*A->cmap->n);CHKERRQ(ierr); 721 PetscFunctionReturn(0); 722 } 723 724 /* -----------------------------------------------------------------*/ 725 #undef __FUNCT__ 726 #define __FUNCT__ "MatGetRow_SeqDense" 727 PetscErrorCode MatGetRow_SeqDense(Mat A,PetscInt row,PetscInt *ncols,PetscInt **cols,PetscScalar **vals) 728 { 729 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 730 PetscScalar *v; 731 PetscErrorCode ierr; 732 PetscInt i; 733 734 PetscFunctionBegin; 735 *ncols = A->cmap->n; 736 if (cols) { 737 ierr = PetscMalloc1(A->cmap->n+1,cols);CHKERRQ(ierr); 738 for (i=0; i<A->cmap->n; i++) (*cols)[i] = i; 739 } 740 if (vals) { 741 ierr = PetscMalloc1(A->cmap->n+1,vals);CHKERRQ(ierr); 742 v = mat->v + row; 743 for (i=0; i<A->cmap->n; i++) {(*vals)[i] = *v; v += mat->lda;} 744 } 745 PetscFunctionReturn(0); 746 } 747 748 #undef __FUNCT__ 749 #define __FUNCT__ "MatRestoreRow_SeqDense" 750 PetscErrorCode MatRestoreRow_SeqDense(Mat A,PetscInt row,PetscInt *ncols,PetscInt **cols,PetscScalar **vals) 751 { 752 PetscErrorCode ierr; 753 754 PetscFunctionBegin; 755 if (cols) {ierr = PetscFree(*cols);CHKERRQ(ierr);} 756 if (vals) {ierr = PetscFree(*vals);CHKERRQ(ierr); } 757 PetscFunctionReturn(0); 758 } 759 /* ----------------------------------------------------------------*/ 760 #undef __FUNCT__ 761 #define __FUNCT__ "MatSetValues_SeqDense" 762 PetscErrorCode MatSetValues_SeqDense(Mat A,PetscInt m,const PetscInt indexm[],PetscInt n,const PetscInt indexn[],const PetscScalar v[],InsertMode addv) 763 { 764 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 765 PetscInt i,j,idx=0; 766 767 PetscFunctionBegin; 768 if (!mat->roworiented) { 769 if (addv == INSERT_VALUES) { 770 for (j=0; j<n; j++) { 771 if (indexn[j] < 0) {idx += m; continue;} 772 #if defined(PETSC_USE_DEBUG) 773 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); 774 #endif 775 for (i=0; i<m; i++) { 776 if (indexm[i] < 0) {idx++; continue;} 777 #if defined(PETSC_USE_DEBUG) 778 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); 779 #endif 780 mat->v[indexn[j]*mat->lda + indexm[i]] = v[idx++]; 781 } 782 } 783 } else { 784 for (j=0; j<n; j++) { 785 if (indexn[j] < 0) {idx += m; continue;} 786 #if defined(PETSC_USE_DEBUG) 787 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); 788 #endif 789 for (i=0; i<m; i++) { 790 if (indexm[i] < 0) {idx++; continue;} 791 #if defined(PETSC_USE_DEBUG) 792 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); 793 #endif 794 mat->v[indexn[j]*mat->lda + indexm[i]] += v[idx++]; 795 } 796 } 797 } 798 } else { 799 if (addv == INSERT_VALUES) { 800 for (i=0; i<m; i++) { 801 if (indexm[i] < 0) { idx += n; continue;} 802 #if defined(PETSC_USE_DEBUG) 803 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); 804 #endif 805 for (j=0; j<n; j++) { 806 if (indexn[j] < 0) { idx++; continue;} 807 #if defined(PETSC_USE_DEBUG) 808 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); 809 #endif 810 mat->v[indexn[j]*mat->lda + indexm[i]] = v[idx++]; 811 } 812 } 813 } else { 814 for (i=0; i<m; i++) { 815 if (indexm[i] < 0) { idx += n; continue;} 816 #if defined(PETSC_USE_DEBUG) 817 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); 818 #endif 819 for (j=0; j<n; j++) { 820 if (indexn[j] < 0) { idx++; continue;} 821 #if defined(PETSC_USE_DEBUG) 822 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); 823 #endif 824 mat->v[indexn[j]*mat->lda + indexm[i]] += v[idx++]; 825 } 826 } 827 } 828 } 829 PetscFunctionReturn(0); 830 } 831 832 #undef __FUNCT__ 833 #define __FUNCT__ "MatGetValues_SeqDense" 834 PetscErrorCode MatGetValues_SeqDense(Mat A,PetscInt m,const PetscInt indexm[],PetscInt n,const PetscInt indexn[],PetscScalar v[]) 835 { 836 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 837 PetscInt i,j; 838 839 PetscFunctionBegin; 840 /* row-oriented output */ 841 for (i=0; i<m; i++) { 842 if (indexm[i] < 0) {v += n;continue;} 843 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); 844 for (j=0; j<n; j++) { 845 if (indexn[j] < 0) {v++; continue;} 846 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); 847 *v++ = mat->v[indexn[j]*mat->lda + indexm[i]]; 848 } 849 } 850 PetscFunctionReturn(0); 851 } 852 853 /* -----------------------------------------------------------------*/ 854 855 #undef __FUNCT__ 856 #define __FUNCT__ "MatLoad_SeqDense" 857 PetscErrorCode MatLoad_SeqDense(Mat newmat,PetscViewer viewer) 858 { 859 Mat_SeqDense *a; 860 PetscErrorCode ierr; 861 PetscInt *scols,i,j,nz,header[4]; 862 int fd; 863 PetscMPIInt size; 864 PetscInt *rowlengths = 0,M,N,*cols,grows,gcols; 865 PetscScalar *vals,*svals,*v,*w; 866 MPI_Comm comm; 867 868 PetscFunctionBegin; 869 /* force binary viewer to load .info file if it has not yet done so */ 870 ierr = PetscViewerSetUp(viewer);CHKERRQ(ierr); 871 ierr = PetscObjectGetComm((PetscObject)viewer,&comm);CHKERRQ(ierr); 872 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 873 if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"view must have one processor"); 874 ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 875 ierr = PetscBinaryRead(fd,header,4,PETSC_INT);CHKERRQ(ierr); 876 if (header[0] != MAT_FILE_CLASSID) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"Not matrix object"); 877 M = header[1]; N = header[2]; nz = header[3]; 878 879 /* set global size if not set already*/ 880 if (newmat->rmap->n < 0 && newmat->rmap->N < 0 && newmat->cmap->n < 0 && newmat->cmap->N < 0) { 881 ierr = MatSetSizes(newmat,M,N,M,N);CHKERRQ(ierr); 882 } else { 883 /* if sizes and type are already set, check if the vector global sizes are correct */ 884 ierr = MatGetSize(newmat,&grows,&gcols);CHKERRQ(ierr); 885 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); 886 } 887 a = (Mat_SeqDense*)newmat->data; 888 if (!a->user_alloc) { 889 ierr = MatSeqDenseSetPreallocation(newmat,NULL);CHKERRQ(ierr); 890 } 891 892 if (nz == MATRIX_BINARY_FORMAT_DENSE) { /* matrix in file is dense */ 893 a = (Mat_SeqDense*)newmat->data; 894 v = a->v; 895 /* Allocate some temp space to read in the values and then flip them 896 from row major to column major */ 897 ierr = PetscMalloc1(M*N > 0 ? M*N : 1,&w);CHKERRQ(ierr); 898 /* read in nonzero values */ 899 ierr = PetscBinaryRead(fd,w,M*N,PETSC_SCALAR);CHKERRQ(ierr); 900 /* now flip the values and store them in the matrix*/ 901 for (j=0; j<N; j++) { 902 for (i=0; i<M; i++) { 903 *v++ =w[i*N+j]; 904 } 905 } 906 ierr = PetscFree(w);CHKERRQ(ierr); 907 } else { 908 /* read row lengths */ 909 ierr = PetscMalloc1(M+1,&rowlengths);CHKERRQ(ierr); 910 ierr = PetscBinaryRead(fd,rowlengths,M,PETSC_INT);CHKERRQ(ierr); 911 912 a = (Mat_SeqDense*)newmat->data; 913 v = a->v; 914 915 /* read column indices and nonzeros */ 916 ierr = PetscMalloc1(nz+1,&scols);CHKERRQ(ierr); 917 cols = scols; 918 ierr = PetscBinaryRead(fd,cols,nz,PETSC_INT);CHKERRQ(ierr); 919 ierr = PetscMalloc1(nz+1,&svals);CHKERRQ(ierr); 920 vals = svals; 921 ierr = PetscBinaryRead(fd,vals,nz,PETSC_SCALAR);CHKERRQ(ierr); 922 923 /* insert into matrix */ 924 for (i=0; i<M; i++) { 925 for (j=0; j<rowlengths[i]; j++) v[i+M*scols[j]] = svals[j]; 926 svals += rowlengths[i]; scols += rowlengths[i]; 927 } 928 ierr = PetscFree(vals);CHKERRQ(ierr); 929 ierr = PetscFree(cols);CHKERRQ(ierr); 930 ierr = PetscFree(rowlengths);CHKERRQ(ierr); 931 } 932 ierr = MatAssemblyBegin(newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 933 ierr = MatAssemblyEnd(newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 934 PetscFunctionReturn(0); 935 } 936 937 #undef __FUNCT__ 938 #define __FUNCT__ "MatView_SeqDense_ASCII" 939 static PetscErrorCode MatView_SeqDense_ASCII(Mat A,PetscViewer viewer) 940 { 941 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 942 PetscErrorCode ierr; 943 PetscInt i,j; 944 const char *name; 945 PetscScalar *v; 946 PetscViewerFormat format; 947 #if defined(PETSC_USE_COMPLEX) 948 PetscBool allreal = PETSC_TRUE; 949 #endif 950 951 PetscFunctionBegin; 952 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 953 if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) { 954 PetscFunctionReturn(0); /* do nothing for now */ 955 } else if (format == PETSC_VIEWER_ASCII_COMMON) { 956 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);CHKERRQ(ierr); 957 for (i=0; i<A->rmap->n; i++) { 958 v = a->v + i; 959 ierr = PetscViewerASCIIPrintf(viewer,"row %D:",i);CHKERRQ(ierr); 960 for (j=0; j<A->cmap->n; j++) { 961 #if defined(PETSC_USE_COMPLEX) 962 if (PetscRealPart(*v) != 0.0 && PetscImaginaryPart(*v) != 0.0) { 963 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g + %g i) ",j,(double)PetscRealPart(*v),(double)PetscImaginaryPart(*v));CHKERRQ(ierr); 964 } else if (PetscRealPart(*v)) { 965 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g) ",j,(double)PetscRealPart(*v));CHKERRQ(ierr); 966 } 967 #else 968 if (*v) { 969 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g) ",j,(double)*v);CHKERRQ(ierr); 970 } 971 #endif 972 v += a->lda; 973 } 974 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 975 } 976 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);CHKERRQ(ierr); 977 } else { 978 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);CHKERRQ(ierr); 979 #if defined(PETSC_USE_COMPLEX) 980 /* determine if matrix has all real values */ 981 v = a->v; 982 for (i=0; i<A->rmap->n*A->cmap->n; i++) { 983 if (PetscImaginaryPart(v[i])) { allreal = PETSC_FALSE; break;} 984 } 985 #endif 986 if (format == PETSC_VIEWER_ASCII_MATLAB) { 987 ierr = PetscObjectGetName((PetscObject)A,&name);CHKERRQ(ierr); 988 ierr = PetscViewerASCIIPrintf(viewer,"%% Size = %D %D \n",A->rmap->n,A->cmap->n);CHKERRQ(ierr); 989 ierr = PetscViewerASCIIPrintf(viewer,"%s = zeros(%D,%D);\n",name,A->rmap->n,A->cmap->n);CHKERRQ(ierr); 990 ierr = PetscViewerASCIIPrintf(viewer,"%s = [\n",name);CHKERRQ(ierr); 991 } 992 993 for (i=0; i<A->rmap->n; i++) { 994 v = a->v + i; 995 for (j=0; j<A->cmap->n; j++) { 996 #if defined(PETSC_USE_COMPLEX) 997 if (allreal) { 998 ierr = PetscViewerASCIIPrintf(viewer,"%18.16e ",(double)PetscRealPart(*v));CHKERRQ(ierr); 999 } else { 1000 ierr = PetscViewerASCIIPrintf(viewer,"%18.16e + %18.16ei ",(double)PetscRealPart(*v),(double)PetscImaginaryPart(*v));CHKERRQ(ierr); 1001 } 1002 #else 1003 ierr = PetscViewerASCIIPrintf(viewer,"%18.16e ",(double)*v);CHKERRQ(ierr); 1004 #endif 1005 v += a->lda; 1006 } 1007 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 1008 } 1009 if (format == PETSC_VIEWER_ASCII_MATLAB) { 1010 ierr = PetscViewerASCIIPrintf(viewer,"];\n");CHKERRQ(ierr); 1011 } 1012 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);CHKERRQ(ierr); 1013 } 1014 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 1015 PetscFunctionReturn(0); 1016 } 1017 1018 #undef __FUNCT__ 1019 #define __FUNCT__ "MatView_SeqDense_Binary" 1020 static PetscErrorCode MatView_SeqDense_Binary(Mat A,PetscViewer viewer) 1021 { 1022 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1023 PetscErrorCode ierr; 1024 int fd; 1025 PetscInt ict,j,n = A->cmap->n,m = A->rmap->n,i,*col_lens,nz = m*n; 1026 PetscScalar *v,*anonz,*vals; 1027 PetscViewerFormat format; 1028 1029 PetscFunctionBegin; 1030 ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 1031 1032 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 1033 if (format == PETSC_VIEWER_NATIVE) { 1034 /* store the matrix as a dense matrix */ 1035 ierr = PetscMalloc1(4,&col_lens);CHKERRQ(ierr); 1036 1037 col_lens[0] = MAT_FILE_CLASSID; 1038 col_lens[1] = m; 1039 col_lens[2] = n; 1040 col_lens[3] = MATRIX_BINARY_FORMAT_DENSE; 1041 1042 ierr = PetscBinaryWrite(fd,col_lens,4,PETSC_INT,PETSC_TRUE);CHKERRQ(ierr); 1043 ierr = PetscFree(col_lens);CHKERRQ(ierr); 1044 1045 /* write out matrix, by rows */ 1046 ierr = PetscMalloc1(m*n+1,&vals);CHKERRQ(ierr); 1047 v = a->v; 1048 for (j=0; j<n; j++) { 1049 for (i=0; i<m; i++) { 1050 vals[j + i*n] = *v++; 1051 } 1052 } 1053 ierr = PetscBinaryWrite(fd,vals,n*m,PETSC_SCALAR,PETSC_FALSE);CHKERRQ(ierr); 1054 ierr = PetscFree(vals);CHKERRQ(ierr); 1055 } else { 1056 ierr = PetscMalloc1(4+nz,&col_lens);CHKERRQ(ierr); 1057 1058 col_lens[0] = MAT_FILE_CLASSID; 1059 col_lens[1] = m; 1060 col_lens[2] = n; 1061 col_lens[3] = nz; 1062 1063 /* store lengths of each row and write (including header) to file */ 1064 for (i=0; i<m; i++) col_lens[4+i] = n; 1065 ierr = PetscBinaryWrite(fd,col_lens,4+m,PETSC_INT,PETSC_TRUE);CHKERRQ(ierr); 1066 1067 /* Possibly should write in smaller increments, not whole matrix at once? */ 1068 /* store column indices (zero start index) */ 1069 ict = 0; 1070 for (i=0; i<m; i++) { 1071 for (j=0; j<n; j++) col_lens[ict++] = j; 1072 } 1073 ierr = PetscBinaryWrite(fd,col_lens,nz,PETSC_INT,PETSC_FALSE);CHKERRQ(ierr); 1074 ierr = PetscFree(col_lens);CHKERRQ(ierr); 1075 1076 /* store nonzero values */ 1077 ierr = PetscMalloc1(nz+1,&anonz);CHKERRQ(ierr); 1078 ict = 0; 1079 for (i=0; i<m; i++) { 1080 v = a->v + i; 1081 for (j=0; j<n; j++) { 1082 anonz[ict++] = *v; v += a->lda; 1083 } 1084 } 1085 ierr = PetscBinaryWrite(fd,anonz,nz,PETSC_SCALAR,PETSC_FALSE);CHKERRQ(ierr); 1086 ierr = PetscFree(anonz);CHKERRQ(ierr); 1087 } 1088 PetscFunctionReturn(0); 1089 } 1090 1091 #include <petscdraw.h> 1092 #undef __FUNCT__ 1093 #define __FUNCT__ "MatView_SeqDense_Draw_Zoom" 1094 PetscErrorCode MatView_SeqDense_Draw_Zoom(PetscDraw draw,void *Aa) 1095 { 1096 Mat A = (Mat) Aa; 1097 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1098 PetscErrorCode ierr; 1099 PetscInt m = A->rmap->n,n = A->cmap->n,color,i,j; 1100 PetscScalar *v = a->v; 1101 PetscViewer viewer; 1102 PetscDraw popup; 1103 PetscReal xl,yl,xr,yr,x_l,x_r,y_l,y_r,scale,maxv = 0.0; 1104 PetscViewerFormat format; 1105 1106 PetscFunctionBegin; 1107 ierr = PetscObjectQuery((PetscObject)A,"Zoomviewer",(PetscObject*)&viewer);CHKERRQ(ierr); 1108 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 1109 ierr = PetscDrawGetCoordinates(draw,&xl,&yl,&xr,&yr);CHKERRQ(ierr); 1110 1111 /* Loop over matrix elements drawing boxes */ 1112 if (format != PETSC_VIEWER_DRAW_CONTOUR) { 1113 /* Blue for negative and Red for positive */ 1114 color = PETSC_DRAW_BLUE; 1115 for (j = 0; j < n; j++) { 1116 x_l = j; 1117 x_r = x_l + 1.0; 1118 for (i = 0; i < m; i++) { 1119 y_l = m - i - 1.0; 1120 y_r = y_l + 1.0; 1121 if (PetscRealPart(v[j*m+i]) > 0.) { 1122 color = PETSC_DRAW_RED; 1123 } else if (PetscRealPart(v[j*m+i]) < 0.) { 1124 color = PETSC_DRAW_BLUE; 1125 } else { 1126 continue; 1127 } 1128 ierr = PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);CHKERRQ(ierr); 1129 } 1130 } 1131 } else { 1132 /* use contour shading to indicate magnitude of values */ 1133 /* first determine max of all nonzero values */ 1134 for (i = 0; i < m*n; i++) { 1135 if (PetscAbsScalar(v[i]) > maxv) maxv = PetscAbsScalar(v[i]); 1136 } 1137 scale = (245.0 - PETSC_DRAW_BASIC_COLORS)/maxv; 1138 ierr = PetscDrawGetPopup(draw,&popup);CHKERRQ(ierr); 1139 if (popup) {ierr = PetscDrawScalePopup(popup,0.0,maxv);CHKERRQ(ierr);} 1140 for (j = 0; j < n; j++) { 1141 x_l = j; 1142 x_r = x_l + 1.0; 1143 for (i = 0; i < m; i++) { 1144 y_l = m - i - 1.0; 1145 y_r = y_l + 1.0; 1146 color = PETSC_DRAW_BASIC_COLORS + (int)(scale*PetscAbsScalar(v[j*m+i])); 1147 ierr = PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);CHKERRQ(ierr); 1148 } 1149 } 1150 } 1151 PetscFunctionReturn(0); 1152 } 1153 1154 #undef __FUNCT__ 1155 #define __FUNCT__ "MatView_SeqDense_Draw" 1156 PetscErrorCode MatView_SeqDense_Draw(Mat A,PetscViewer viewer) 1157 { 1158 PetscDraw draw; 1159 PetscBool isnull; 1160 PetscReal xr,yr,xl,yl,h,w; 1161 PetscErrorCode ierr; 1162 1163 PetscFunctionBegin; 1164 ierr = PetscViewerDrawGetDraw(viewer,0,&draw);CHKERRQ(ierr); 1165 ierr = PetscDrawIsNull(draw,&isnull);CHKERRQ(ierr); 1166 if (isnull) PetscFunctionReturn(0); 1167 1168 ierr = PetscObjectCompose((PetscObject)A,"Zoomviewer",(PetscObject)viewer);CHKERRQ(ierr); 1169 xr = A->cmap->n; yr = A->rmap->n; h = yr/10.0; w = xr/10.0; 1170 xr += w; yr += h; xl = -w; yl = -h; 1171 ierr = PetscDrawSetCoordinates(draw,xl,yl,xr,yr);CHKERRQ(ierr); 1172 ierr = PetscDrawZoom(draw,MatView_SeqDense_Draw_Zoom,A);CHKERRQ(ierr); 1173 ierr = PetscObjectCompose((PetscObject)A,"Zoomviewer",NULL);CHKERRQ(ierr); 1174 PetscFunctionReturn(0); 1175 } 1176 1177 #undef __FUNCT__ 1178 #define __FUNCT__ "MatView_SeqDense" 1179 PetscErrorCode MatView_SeqDense(Mat A,PetscViewer viewer) 1180 { 1181 PetscErrorCode ierr; 1182 PetscBool iascii,isbinary,isdraw; 1183 1184 PetscFunctionBegin; 1185 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 1186 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);CHKERRQ(ierr); 1187 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);CHKERRQ(ierr); 1188 1189 if (iascii) { 1190 ierr = MatView_SeqDense_ASCII(A,viewer);CHKERRQ(ierr); 1191 } else if (isbinary) { 1192 ierr = MatView_SeqDense_Binary(A,viewer);CHKERRQ(ierr); 1193 } else if (isdraw) { 1194 ierr = MatView_SeqDense_Draw(A,viewer);CHKERRQ(ierr); 1195 } 1196 PetscFunctionReturn(0); 1197 } 1198 1199 #undef __FUNCT__ 1200 #define __FUNCT__ "MatDestroy_SeqDense" 1201 PetscErrorCode MatDestroy_SeqDense(Mat mat) 1202 { 1203 Mat_SeqDense *l = (Mat_SeqDense*)mat->data; 1204 PetscErrorCode ierr; 1205 1206 PetscFunctionBegin; 1207 #if defined(PETSC_USE_LOG) 1208 PetscLogObjectState((PetscObject)mat,"Rows %D Cols %D",mat->rmap->n,mat->cmap->n); 1209 #endif 1210 ierr = PetscFree(l->pivots);CHKERRQ(ierr); 1211 if (!l->user_alloc) {ierr = PetscFree(l->v);CHKERRQ(ierr);} 1212 ierr = PetscFree(mat->data);CHKERRQ(ierr); 1213 1214 ierr = PetscObjectChangeTypeName((PetscObject)mat,0);CHKERRQ(ierr); 1215 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetArray_C",NULL);CHKERRQ(ierr); 1216 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatDenseRestoreArray_C",NULL);CHKERRQ(ierr); 1217 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatConvert_seqdense_seqaij_C",NULL);CHKERRQ(ierr); 1218 #if defined(PETSC_HAVE_ELEMENTAL) 1219 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatConvert_seqdense_elemental_C",NULL);CHKERRQ(ierr); 1220 #endif 1221 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatSeqDenseSetPreallocation_C",NULL);CHKERRQ(ierr); 1222 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatMatMult_seqaij_seqdense_C",NULL);CHKERRQ(ierr); 1223 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatMatMultSymbolic_seqaij_seqdense_C",NULL);CHKERRQ(ierr); 1224 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatMatMultNumeric_seqaij_seqdense_C",NULL);CHKERRQ(ierr); 1225 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatTransposeMatMult_seqaij_seqdense_C",NULL);CHKERRQ(ierr); 1226 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatTransposeMatMultSymbolic_seqaij_seqdense_C",NULL);CHKERRQ(ierr); 1227 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatTransposeMatMultNumeric_seqaij_seqdense_C",NULL);CHKERRQ(ierr); 1228 PetscFunctionReturn(0); 1229 } 1230 1231 #undef __FUNCT__ 1232 #define __FUNCT__ "MatTranspose_SeqDense" 1233 PetscErrorCode MatTranspose_SeqDense(Mat A,MatReuse reuse,Mat *matout) 1234 { 1235 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1236 PetscErrorCode ierr; 1237 PetscInt k,j,m,n,M; 1238 PetscScalar *v,tmp; 1239 1240 PetscFunctionBegin; 1241 v = mat->v; m = A->rmap->n; M = mat->lda; n = A->cmap->n; 1242 if (reuse == MAT_REUSE_MATRIX && *matout == A) { /* in place transpose */ 1243 if (m != n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Can not transpose non-square matrix in place"); 1244 else { 1245 for (j=0; j<m; j++) { 1246 for (k=0; k<j; k++) { 1247 tmp = v[j + k*M]; 1248 v[j + k*M] = v[k + j*M]; 1249 v[k + j*M] = tmp; 1250 } 1251 } 1252 } 1253 } else { /* out-of-place transpose */ 1254 Mat tmat; 1255 Mat_SeqDense *tmatd; 1256 PetscScalar *v2; 1257 PetscInt M2; 1258 1259 if (reuse == MAT_INITIAL_MATRIX) { 1260 ierr = MatCreate(PetscObjectComm((PetscObject)A),&tmat);CHKERRQ(ierr); 1261 ierr = MatSetSizes(tmat,A->cmap->n,A->rmap->n,A->cmap->n,A->rmap->n);CHKERRQ(ierr); 1262 ierr = MatSetType(tmat,((PetscObject)A)->type_name);CHKERRQ(ierr); 1263 ierr = MatSeqDenseSetPreallocation(tmat,NULL);CHKERRQ(ierr); 1264 } else { 1265 tmat = *matout; 1266 } 1267 tmatd = (Mat_SeqDense*)tmat->data; 1268 v = mat->v; v2 = tmatd->v; M2 = tmatd->lda; 1269 for (j=0; j<n; j++) { 1270 for (k=0; k<m; k++) v2[j + k*M2] = v[k + j*M]; 1271 } 1272 ierr = MatAssemblyBegin(tmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1273 ierr = MatAssemblyEnd(tmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1274 *matout = tmat; 1275 } 1276 PetscFunctionReturn(0); 1277 } 1278 1279 #undef __FUNCT__ 1280 #define __FUNCT__ "MatEqual_SeqDense" 1281 PetscErrorCode MatEqual_SeqDense(Mat A1,Mat A2,PetscBool *flg) 1282 { 1283 Mat_SeqDense *mat1 = (Mat_SeqDense*)A1->data; 1284 Mat_SeqDense *mat2 = (Mat_SeqDense*)A2->data; 1285 PetscInt i,j; 1286 PetscScalar *v1,*v2; 1287 1288 PetscFunctionBegin; 1289 if (A1->rmap->n != A2->rmap->n) {*flg = PETSC_FALSE; PetscFunctionReturn(0);} 1290 if (A1->cmap->n != A2->cmap->n) {*flg = PETSC_FALSE; PetscFunctionReturn(0);} 1291 for (i=0; i<A1->rmap->n; i++) { 1292 v1 = mat1->v+i; v2 = mat2->v+i; 1293 for (j=0; j<A1->cmap->n; j++) { 1294 if (*v1 != *v2) {*flg = PETSC_FALSE; PetscFunctionReturn(0);} 1295 v1 += mat1->lda; v2 += mat2->lda; 1296 } 1297 } 1298 *flg = PETSC_TRUE; 1299 PetscFunctionReturn(0); 1300 } 1301 1302 #undef __FUNCT__ 1303 #define __FUNCT__ "MatGetDiagonal_SeqDense" 1304 PetscErrorCode MatGetDiagonal_SeqDense(Mat A,Vec v) 1305 { 1306 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1307 PetscErrorCode ierr; 1308 PetscInt i,n,len; 1309 PetscScalar *x,zero = 0.0; 1310 1311 PetscFunctionBegin; 1312 ierr = VecSet(v,zero);CHKERRQ(ierr); 1313 ierr = VecGetSize(v,&n);CHKERRQ(ierr); 1314 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 1315 len = PetscMin(A->rmap->n,A->cmap->n); 1316 if (n != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming mat and vec"); 1317 for (i=0; i<len; i++) { 1318 x[i] = mat->v[i*mat->lda + i]; 1319 } 1320 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 1321 PetscFunctionReturn(0); 1322 } 1323 1324 #undef __FUNCT__ 1325 #define __FUNCT__ "MatDiagonalScale_SeqDense" 1326 PetscErrorCode MatDiagonalScale_SeqDense(Mat A,Vec ll,Vec rr) 1327 { 1328 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1329 const PetscScalar *l,*r; 1330 PetscScalar x,*v; 1331 PetscErrorCode ierr; 1332 PetscInt i,j,m = A->rmap->n,n = A->cmap->n; 1333 1334 PetscFunctionBegin; 1335 if (ll) { 1336 ierr = VecGetSize(ll,&m);CHKERRQ(ierr); 1337 ierr = VecGetArrayRead(ll,&l);CHKERRQ(ierr); 1338 if (m != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Left scaling vec wrong size"); 1339 for (i=0; i<m; i++) { 1340 x = l[i]; 1341 v = mat->v + i; 1342 for (j=0; j<n; j++) { (*v) *= x; v+= m;} 1343 } 1344 ierr = VecRestoreArrayRead(ll,&l);CHKERRQ(ierr); 1345 ierr = PetscLogFlops(1.0*n*m);CHKERRQ(ierr); 1346 } 1347 if (rr) { 1348 ierr = VecGetSize(rr,&n);CHKERRQ(ierr); 1349 ierr = VecGetArrayRead(rr,&r);CHKERRQ(ierr); 1350 if (n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Right scaling vec wrong size"); 1351 for (i=0; i<n; i++) { 1352 x = r[i]; 1353 v = mat->v + i*m; 1354 for (j=0; j<m; j++) (*v++) *= x; 1355 } 1356 ierr = VecRestoreArrayRead(rr,&r);CHKERRQ(ierr); 1357 ierr = PetscLogFlops(1.0*n*m);CHKERRQ(ierr); 1358 } 1359 PetscFunctionReturn(0); 1360 } 1361 1362 #undef __FUNCT__ 1363 #define __FUNCT__ "MatNorm_SeqDense" 1364 PetscErrorCode MatNorm_SeqDense(Mat A,NormType type,PetscReal *nrm) 1365 { 1366 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1367 PetscScalar *v = mat->v; 1368 PetscReal sum = 0.0; 1369 PetscInt lda =mat->lda,m=A->rmap->n,i,j; 1370 PetscErrorCode ierr; 1371 1372 PetscFunctionBegin; 1373 if (type == NORM_FROBENIUS) { 1374 if (lda>m) { 1375 for (j=0; j<A->cmap->n; j++) { 1376 v = mat->v+j*lda; 1377 for (i=0; i<m; i++) { 1378 sum += PetscRealPart(PetscConj(*v)*(*v)); v++; 1379 } 1380 } 1381 } else { 1382 for (i=0; i<A->cmap->n*A->rmap->n; i++) { 1383 sum += PetscRealPart(PetscConj(*v)*(*v)); v++; 1384 } 1385 } 1386 *nrm = PetscSqrtReal(sum); 1387 ierr = PetscLogFlops(2.0*A->cmap->n*A->rmap->n);CHKERRQ(ierr); 1388 } else if (type == NORM_1) { 1389 *nrm = 0.0; 1390 for (j=0; j<A->cmap->n; j++) { 1391 v = mat->v + j*mat->lda; 1392 sum = 0.0; 1393 for (i=0; i<A->rmap->n; i++) { 1394 sum += PetscAbsScalar(*v); v++; 1395 } 1396 if (sum > *nrm) *nrm = sum; 1397 } 1398 ierr = PetscLogFlops(1.0*A->cmap->n*A->rmap->n);CHKERRQ(ierr); 1399 } else if (type == NORM_INFINITY) { 1400 *nrm = 0.0; 1401 for (j=0; j<A->rmap->n; j++) { 1402 v = mat->v + j; 1403 sum = 0.0; 1404 for (i=0; i<A->cmap->n; i++) { 1405 sum += PetscAbsScalar(*v); v += mat->lda; 1406 } 1407 if (sum > *nrm) *nrm = sum; 1408 } 1409 ierr = PetscLogFlops(1.0*A->cmap->n*A->rmap->n);CHKERRQ(ierr); 1410 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No two norm"); 1411 PetscFunctionReturn(0); 1412 } 1413 1414 #undef __FUNCT__ 1415 #define __FUNCT__ "MatSetOption_SeqDense" 1416 PetscErrorCode MatSetOption_SeqDense(Mat A,MatOption op,PetscBool flg) 1417 { 1418 Mat_SeqDense *aij = (Mat_SeqDense*)A->data; 1419 PetscErrorCode ierr; 1420 1421 PetscFunctionBegin; 1422 switch (op) { 1423 case MAT_ROW_ORIENTED: 1424 aij->roworiented = flg; 1425 break; 1426 case MAT_NEW_NONZERO_LOCATIONS: 1427 case MAT_NEW_NONZERO_LOCATION_ERR: 1428 case MAT_NEW_NONZERO_ALLOCATION_ERR: 1429 case MAT_NEW_DIAGONALS: 1430 case MAT_KEEP_NONZERO_PATTERN: 1431 case MAT_IGNORE_OFF_PROC_ENTRIES: 1432 case MAT_USE_HASH_TABLE: 1433 case MAT_IGNORE_LOWER_TRIANGULAR: 1434 ierr = PetscInfo1(A,"Option %s ignored\n",MatOptions[op]);CHKERRQ(ierr); 1435 break; 1436 case MAT_SPD: 1437 case MAT_SYMMETRIC: 1438 case MAT_STRUCTURALLY_SYMMETRIC: 1439 case MAT_HERMITIAN: 1440 case MAT_SYMMETRY_ETERNAL: 1441 /* These options are handled directly by MatSetOption() */ 1442 break; 1443 default: 1444 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"unknown option %s",MatOptions[op]); 1445 } 1446 PetscFunctionReturn(0); 1447 } 1448 1449 #undef __FUNCT__ 1450 #define __FUNCT__ "MatZeroEntries_SeqDense" 1451 PetscErrorCode MatZeroEntries_SeqDense(Mat A) 1452 { 1453 Mat_SeqDense *l = (Mat_SeqDense*)A->data; 1454 PetscErrorCode ierr; 1455 PetscInt lda=l->lda,m=A->rmap->n,j; 1456 1457 PetscFunctionBegin; 1458 if (lda>m) { 1459 for (j=0; j<A->cmap->n; j++) { 1460 ierr = PetscMemzero(l->v+j*lda,m*sizeof(PetscScalar));CHKERRQ(ierr); 1461 } 1462 } else { 1463 ierr = PetscMemzero(l->v,A->rmap->n*A->cmap->n*sizeof(PetscScalar));CHKERRQ(ierr); 1464 } 1465 PetscFunctionReturn(0); 1466 } 1467 1468 #undef __FUNCT__ 1469 #define __FUNCT__ "MatZeroRows_SeqDense" 1470 PetscErrorCode MatZeroRows_SeqDense(Mat A,PetscInt N,const PetscInt rows[],PetscScalar diag,Vec x,Vec b) 1471 { 1472 PetscErrorCode ierr; 1473 Mat_SeqDense *l = (Mat_SeqDense*)A->data; 1474 PetscInt m = l->lda, n = A->cmap->n, i,j; 1475 PetscScalar *slot,*bb; 1476 const PetscScalar *xx; 1477 1478 PetscFunctionBegin; 1479 #if defined(PETSC_USE_DEBUG) 1480 for (i=0; i<N; i++) { 1481 if (rows[i] < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative row requested to be zeroed"); 1482 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); 1483 } 1484 #endif 1485 1486 /* fix right hand side if needed */ 1487 if (x && b) { 1488 ierr = VecGetArrayRead(x,&xx);CHKERRQ(ierr); 1489 ierr = VecGetArray(b,&bb);CHKERRQ(ierr); 1490 for (i=0; i<N; i++) bb[rows[i]] = diag*xx[rows[i]]; 1491 ierr = VecRestoreArrayRead(x,&xx);CHKERRQ(ierr); 1492 ierr = VecRestoreArray(b,&bb);CHKERRQ(ierr); 1493 } 1494 1495 for (i=0; i<N; i++) { 1496 slot = l->v + rows[i]; 1497 for (j=0; j<n; j++) { *slot = 0.0; slot += m;} 1498 } 1499 if (diag != 0.0) { 1500 if (A->rmap->n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only coded for square matrices"); 1501 for (i=0; i<N; i++) { 1502 slot = l->v + (m+1)*rows[i]; 1503 *slot = diag; 1504 } 1505 } 1506 PetscFunctionReturn(0); 1507 } 1508 1509 #undef __FUNCT__ 1510 #define __FUNCT__ "MatDenseGetArray_SeqDense" 1511 PetscErrorCode MatDenseGetArray_SeqDense(Mat A,PetscScalar *array[]) 1512 { 1513 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1514 1515 PetscFunctionBegin; 1516 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"); 1517 *array = mat->v; 1518 PetscFunctionReturn(0); 1519 } 1520 1521 #undef __FUNCT__ 1522 #define __FUNCT__ "MatDenseRestoreArray_SeqDense" 1523 PetscErrorCode MatDenseRestoreArray_SeqDense(Mat A,PetscScalar *array[]) 1524 { 1525 PetscFunctionBegin; 1526 *array = 0; /* user cannot accidently use the array later */ 1527 PetscFunctionReturn(0); 1528 } 1529 1530 #undef __FUNCT__ 1531 #define __FUNCT__ "MatDenseGetArray" 1532 /*@C 1533 MatDenseGetArray - gives access to the array where the data for a SeqDense matrix is stored 1534 1535 Not Collective 1536 1537 Input Parameter: 1538 . mat - a MATSEQDENSE or MATMPIDENSE matrix 1539 1540 Output Parameter: 1541 . array - pointer to the data 1542 1543 Level: intermediate 1544 1545 .seealso: MatDenseRestoreArray() 1546 @*/ 1547 PetscErrorCode MatDenseGetArray(Mat A,PetscScalar **array) 1548 { 1549 PetscErrorCode ierr; 1550 1551 PetscFunctionBegin; 1552 ierr = PetscUseMethod(A,"MatDenseGetArray_C",(Mat,PetscScalar**),(A,array));CHKERRQ(ierr); 1553 PetscFunctionReturn(0); 1554 } 1555 1556 #undef __FUNCT__ 1557 #define __FUNCT__ "MatDenseRestoreArray" 1558 /*@C 1559 MatDenseRestoreArray - returns access to the array where the data for a dense matrix is stored obtained by MatDenseGetArray() 1560 1561 Not Collective 1562 1563 Input Parameters: 1564 . mat - a MATSEQDENSE or MATMPIDENSE matrix 1565 . array - pointer to the data 1566 1567 Level: intermediate 1568 1569 .seealso: MatDenseGetArray() 1570 @*/ 1571 PetscErrorCode MatDenseRestoreArray(Mat A,PetscScalar **array) 1572 { 1573 PetscErrorCode ierr; 1574 1575 PetscFunctionBegin; 1576 ierr = PetscUseMethod(A,"MatDenseRestoreArray_C",(Mat,PetscScalar**),(A,array));CHKERRQ(ierr); 1577 PetscFunctionReturn(0); 1578 } 1579 1580 #undef __FUNCT__ 1581 #define __FUNCT__ "MatGetSubMatrix_SeqDense" 1582 static PetscErrorCode MatGetSubMatrix_SeqDense(Mat A,IS isrow,IS iscol,PetscInt cs,MatReuse scall,Mat *B) 1583 { 1584 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1585 PetscErrorCode ierr; 1586 PetscInt i,j,nrows,ncols; 1587 const PetscInt *irow,*icol; 1588 PetscScalar *av,*bv,*v = mat->v; 1589 Mat newmat; 1590 1591 PetscFunctionBegin; 1592 ierr = ISGetIndices(isrow,&irow);CHKERRQ(ierr); 1593 ierr = ISGetIndices(iscol,&icol);CHKERRQ(ierr); 1594 ierr = ISGetLocalSize(isrow,&nrows);CHKERRQ(ierr); 1595 ierr = ISGetLocalSize(iscol,&ncols);CHKERRQ(ierr); 1596 1597 /* Check submatrixcall */ 1598 if (scall == MAT_REUSE_MATRIX) { 1599 PetscInt n_cols,n_rows; 1600 ierr = MatGetSize(*B,&n_rows,&n_cols);CHKERRQ(ierr); 1601 if (n_rows != nrows || n_cols != ncols) { 1602 /* resize the result matrix to match number of requested rows/columns */ 1603 ierr = MatSetSizes(*B,nrows,ncols,nrows,ncols);CHKERRQ(ierr); 1604 } 1605 newmat = *B; 1606 } else { 1607 /* Create and fill new matrix */ 1608 ierr = MatCreate(PetscObjectComm((PetscObject)A),&newmat);CHKERRQ(ierr); 1609 ierr = MatSetSizes(newmat,nrows,ncols,nrows,ncols);CHKERRQ(ierr); 1610 ierr = MatSetType(newmat,((PetscObject)A)->type_name);CHKERRQ(ierr); 1611 ierr = MatSeqDenseSetPreallocation(newmat,NULL);CHKERRQ(ierr); 1612 } 1613 1614 /* Now extract the data pointers and do the copy,column at a time */ 1615 bv = ((Mat_SeqDense*)newmat->data)->v; 1616 1617 for (i=0; i<ncols; i++) { 1618 av = v + mat->lda*icol[i]; 1619 for (j=0; j<nrows; j++) *bv++ = av[irow[j]]; 1620 } 1621 1622 /* Assemble the matrices so that the correct flags are set */ 1623 ierr = MatAssemblyBegin(newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1624 ierr = MatAssemblyEnd(newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1625 1626 /* Free work space */ 1627 ierr = ISRestoreIndices(isrow,&irow);CHKERRQ(ierr); 1628 ierr = ISRestoreIndices(iscol,&icol);CHKERRQ(ierr); 1629 *B = newmat; 1630 PetscFunctionReturn(0); 1631 } 1632 1633 #undef __FUNCT__ 1634 #define __FUNCT__ "MatGetSubMatrices_SeqDense" 1635 PetscErrorCode MatGetSubMatrices_SeqDense(Mat A,PetscInt n,const IS irow[],const IS icol[],MatReuse scall,Mat *B[]) 1636 { 1637 PetscErrorCode ierr; 1638 PetscInt i; 1639 1640 PetscFunctionBegin; 1641 if (scall == MAT_INITIAL_MATRIX) { 1642 ierr = PetscMalloc1(n+1,B);CHKERRQ(ierr); 1643 } 1644 1645 for (i=0; i<n; i++) { 1646 ierr = MatGetSubMatrix_SeqDense(A,irow[i],icol[i],PETSC_DECIDE,scall,&(*B)[i]);CHKERRQ(ierr); 1647 } 1648 PetscFunctionReturn(0); 1649 } 1650 1651 #undef __FUNCT__ 1652 #define __FUNCT__ "MatAssemblyBegin_SeqDense" 1653 PetscErrorCode MatAssemblyBegin_SeqDense(Mat mat,MatAssemblyType mode) 1654 { 1655 PetscFunctionBegin; 1656 PetscFunctionReturn(0); 1657 } 1658 1659 #undef __FUNCT__ 1660 #define __FUNCT__ "MatAssemblyEnd_SeqDense" 1661 PetscErrorCode MatAssemblyEnd_SeqDense(Mat mat,MatAssemblyType mode) 1662 { 1663 PetscFunctionBegin; 1664 PetscFunctionReturn(0); 1665 } 1666 1667 #undef __FUNCT__ 1668 #define __FUNCT__ "MatCopy_SeqDense" 1669 PetscErrorCode MatCopy_SeqDense(Mat A,Mat B,MatStructure str) 1670 { 1671 Mat_SeqDense *a = (Mat_SeqDense*)A->data,*b = (Mat_SeqDense*)B->data; 1672 PetscErrorCode ierr; 1673 PetscInt lda1=a->lda,lda2=b->lda, m=A->rmap->n,n=A->cmap->n, j; 1674 1675 PetscFunctionBegin; 1676 /* If the two matrices don't have the same copy implementation, they aren't compatible for fast copy. */ 1677 if (A->ops->copy != B->ops->copy) { 1678 ierr = MatCopy_Basic(A,B,str);CHKERRQ(ierr); 1679 PetscFunctionReturn(0); 1680 } 1681 if (m != B->rmap->n || n != B->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"size(B) != size(A)"); 1682 if (lda1>m || lda2>m) { 1683 for (j=0; j<n; j++) { 1684 ierr = PetscMemcpy(b->v+j*lda2,a->v+j*lda1,m*sizeof(PetscScalar));CHKERRQ(ierr); 1685 } 1686 } else { 1687 ierr = PetscMemcpy(b->v,a->v,A->rmap->n*A->cmap->n*sizeof(PetscScalar));CHKERRQ(ierr); 1688 } 1689 PetscFunctionReturn(0); 1690 } 1691 1692 #undef __FUNCT__ 1693 #define __FUNCT__ "MatSetUp_SeqDense" 1694 PetscErrorCode MatSetUp_SeqDense(Mat A) 1695 { 1696 PetscErrorCode ierr; 1697 1698 PetscFunctionBegin; 1699 ierr = MatSeqDenseSetPreallocation(A,0);CHKERRQ(ierr); 1700 PetscFunctionReturn(0); 1701 } 1702 1703 #undef __FUNCT__ 1704 #define __FUNCT__ "MatConjugate_SeqDense" 1705 static PetscErrorCode MatConjugate_SeqDense(Mat A) 1706 { 1707 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1708 PetscInt i,nz = A->rmap->n*A->cmap->n; 1709 PetscScalar *aa = a->v; 1710 1711 PetscFunctionBegin; 1712 for (i=0; i<nz; i++) aa[i] = PetscConj(aa[i]); 1713 PetscFunctionReturn(0); 1714 } 1715 1716 #undef __FUNCT__ 1717 #define __FUNCT__ "MatRealPart_SeqDense" 1718 static PetscErrorCode MatRealPart_SeqDense(Mat A) 1719 { 1720 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1721 PetscInt i,nz = A->rmap->n*A->cmap->n; 1722 PetscScalar *aa = a->v; 1723 1724 PetscFunctionBegin; 1725 for (i=0; i<nz; i++) aa[i] = PetscRealPart(aa[i]); 1726 PetscFunctionReturn(0); 1727 } 1728 1729 #undef __FUNCT__ 1730 #define __FUNCT__ "MatImaginaryPart_SeqDense" 1731 static PetscErrorCode MatImaginaryPart_SeqDense(Mat A) 1732 { 1733 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1734 PetscInt i,nz = A->rmap->n*A->cmap->n; 1735 PetscScalar *aa = a->v; 1736 1737 PetscFunctionBegin; 1738 for (i=0; i<nz; i++) aa[i] = PetscImaginaryPart(aa[i]); 1739 PetscFunctionReturn(0); 1740 } 1741 1742 /* ----------------------------------------------------------------*/ 1743 #undef __FUNCT__ 1744 #define __FUNCT__ "MatMatMult_SeqDense_SeqDense" 1745 PetscErrorCode MatMatMult_SeqDense_SeqDense(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C) 1746 { 1747 PetscErrorCode ierr; 1748 1749 PetscFunctionBegin; 1750 if (scall == MAT_INITIAL_MATRIX) { 1751 ierr = PetscLogEventBegin(MAT_MatMultSymbolic,A,B,0,0);CHKERRQ(ierr); 1752 ierr = MatMatMultSymbolic_SeqDense_SeqDense(A,B,fill,C);CHKERRQ(ierr); 1753 ierr = PetscLogEventEnd(MAT_MatMultSymbolic,A,B,0,0);CHKERRQ(ierr); 1754 } 1755 ierr = PetscLogEventBegin(MAT_MatMultNumeric,A,B,0,0);CHKERRQ(ierr); 1756 ierr = MatMatMultNumeric_SeqDense_SeqDense(A,B,*C);CHKERRQ(ierr); 1757 ierr = PetscLogEventEnd(MAT_MatMultNumeric,A,B,0,0);CHKERRQ(ierr); 1758 PetscFunctionReturn(0); 1759 } 1760 1761 #undef __FUNCT__ 1762 #define __FUNCT__ "MatMatMultSymbolic_SeqDense_SeqDense" 1763 PetscErrorCode MatMatMultSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat *C) 1764 { 1765 PetscErrorCode ierr; 1766 PetscInt m=A->rmap->n,n=B->cmap->n; 1767 Mat Cmat; 1768 1769 PetscFunctionBegin; 1770 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); 1771 ierr = MatCreate(PETSC_COMM_SELF,&Cmat);CHKERRQ(ierr); 1772 ierr = MatSetSizes(Cmat,m,n,m,n);CHKERRQ(ierr); 1773 ierr = MatSetType(Cmat,MATSEQDENSE);CHKERRQ(ierr); 1774 ierr = MatSeqDenseSetPreallocation(Cmat,NULL);CHKERRQ(ierr); 1775 1776 *C = Cmat; 1777 PetscFunctionReturn(0); 1778 } 1779 1780 #undef __FUNCT__ 1781 #define __FUNCT__ "MatMatMultNumeric_SeqDense_SeqDense" 1782 PetscErrorCode MatMatMultNumeric_SeqDense_SeqDense(Mat A,Mat B,Mat C) 1783 { 1784 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1785 Mat_SeqDense *b = (Mat_SeqDense*)B->data; 1786 Mat_SeqDense *c = (Mat_SeqDense*)C->data; 1787 PetscBLASInt m,n,k; 1788 PetscScalar _DOne=1.0,_DZero=0.0; 1789 PetscErrorCode ierr; 1790 1791 PetscFunctionBegin; 1792 ierr = PetscBLASIntCast(A->rmap->n,&m);CHKERRQ(ierr); 1793 ierr = PetscBLASIntCast(B->cmap->n,&n);CHKERRQ(ierr); 1794 ierr = PetscBLASIntCast(A->cmap->n,&k);CHKERRQ(ierr); 1795 PetscStackCallBLAS("BLASgemv",BLASgemm_("N","N",&m,&n,&k,&_DOne,a->v,&a->lda,b->v,&b->lda,&_DZero,c->v,&c->lda)); 1796 PetscFunctionReturn(0); 1797 } 1798 1799 #undef __FUNCT__ 1800 #define __FUNCT__ "MatTransposeMatMult_SeqDense_SeqDense" 1801 PetscErrorCode MatTransposeMatMult_SeqDense_SeqDense(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C) 1802 { 1803 PetscErrorCode ierr; 1804 1805 PetscFunctionBegin; 1806 if (scall == MAT_INITIAL_MATRIX) { 1807 ierr = PetscLogEventBegin(MAT_TransposeMatMultSymbolic,A,B,0,0);CHKERRQ(ierr); 1808 ierr = MatTransposeMatMultSymbolic_SeqDense_SeqDense(A,B,fill,C);CHKERRQ(ierr); 1809 ierr = PetscLogEventEnd(MAT_TransposeMatMultSymbolic,A,B,0,0);CHKERRQ(ierr); 1810 } 1811 ierr = PetscLogEventBegin(MAT_TransposeMatMultNumeric,A,B,0,0);CHKERRQ(ierr); 1812 ierr = MatTransposeMatMultNumeric_SeqDense_SeqDense(A,B,*C);CHKERRQ(ierr); 1813 ierr = PetscLogEventEnd(MAT_TransposeMatMultNumeric,A,B,0,0);CHKERRQ(ierr); 1814 PetscFunctionReturn(0); 1815 } 1816 1817 #undef __FUNCT__ 1818 #define __FUNCT__ "MatTransposeMatMultSymbolic_SeqDense_SeqDense" 1819 PetscErrorCode MatTransposeMatMultSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat *C) 1820 { 1821 PetscErrorCode ierr; 1822 PetscInt m=A->cmap->n,n=B->cmap->n; 1823 Mat Cmat; 1824 1825 PetscFunctionBegin; 1826 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); 1827 ierr = MatCreate(PETSC_COMM_SELF,&Cmat);CHKERRQ(ierr); 1828 ierr = MatSetSizes(Cmat,m,n,m,n);CHKERRQ(ierr); 1829 ierr = MatSetType(Cmat,MATSEQDENSE);CHKERRQ(ierr); 1830 ierr = MatSeqDenseSetPreallocation(Cmat,NULL);CHKERRQ(ierr); 1831 1832 Cmat->assembled = PETSC_TRUE; 1833 1834 *C = Cmat; 1835 PetscFunctionReturn(0); 1836 } 1837 1838 #undef __FUNCT__ 1839 #define __FUNCT__ "MatTransposeMatMultNumeric_SeqDense_SeqDense" 1840 PetscErrorCode MatTransposeMatMultNumeric_SeqDense_SeqDense(Mat A,Mat B,Mat C) 1841 { 1842 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1843 Mat_SeqDense *b = (Mat_SeqDense*)B->data; 1844 Mat_SeqDense *c = (Mat_SeqDense*)C->data; 1845 PetscBLASInt m,n,k; 1846 PetscScalar _DOne=1.0,_DZero=0.0; 1847 PetscErrorCode ierr; 1848 1849 PetscFunctionBegin; 1850 ierr = PetscBLASIntCast(A->cmap->n,&m);CHKERRQ(ierr); 1851 ierr = PetscBLASIntCast(B->cmap->n,&n);CHKERRQ(ierr); 1852 ierr = PetscBLASIntCast(A->rmap->n,&k);CHKERRQ(ierr); 1853 /* 1854 Note the m and n arguments below are the number rows and columns of A', not A! 1855 */ 1856 PetscStackCallBLAS("BLASgemv",BLASgemm_("T","N",&m,&n,&k,&_DOne,a->v,&a->lda,b->v,&b->lda,&_DZero,c->v,&c->lda)); 1857 PetscFunctionReturn(0); 1858 } 1859 1860 #undef __FUNCT__ 1861 #define __FUNCT__ "MatGetRowMax_SeqDense" 1862 PetscErrorCode MatGetRowMax_SeqDense(Mat A,Vec v,PetscInt idx[]) 1863 { 1864 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1865 PetscErrorCode ierr; 1866 PetscInt i,j,m = A->rmap->n,n = A->cmap->n,p; 1867 PetscScalar *x; 1868 MatScalar *aa = a->v; 1869 1870 PetscFunctionBegin; 1871 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 1872 1873 ierr = VecSet(v,0.0);CHKERRQ(ierr); 1874 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 1875 ierr = VecGetLocalSize(v,&p);CHKERRQ(ierr); 1876 if (p != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 1877 for (i=0; i<m; i++) { 1878 x[i] = aa[i]; if (idx) idx[i] = 0; 1879 for (j=1; j<n; j++) { 1880 if (PetscRealPart(x[i]) < PetscRealPart(aa[i+m*j])) {x[i] = aa[i + m*j]; if (idx) idx[i] = j;} 1881 } 1882 } 1883 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 1884 PetscFunctionReturn(0); 1885 } 1886 1887 #undef __FUNCT__ 1888 #define __FUNCT__ "MatGetRowMaxAbs_SeqDense" 1889 PetscErrorCode MatGetRowMaxAbs_SeqDense(Mat A,Vec v,PetscInt idx[]) 1890 { 1891 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1892 PetscErrorCode ierr; 1893 PetscInt i,j,m = A->rmap->n,n = A->cmap->n,p; 1894 PetscScalar *x; 1895 PetscReal atmp; 1896 MatScalar *aa = a->v; 1897 1898 PetscFunctionBegin; 1899 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 1900 1901 ierr = VecSet(v,0.0);CHKERRQ(ierr); 1902 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 1903 ierr = VecGetLocalSize(v,&p);CHKERRQ(ierr); 1904 if (p != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 1905 for (i=0; i<m; i++) { 1906 x[i] = PetscAbsScalar(aa[i]); 1907 for (j=1; j<n; j++) { 1908 atmp = PetscAbsScalar(aa[i+m*j]); 1909 if (PetscAbsScalar(x[i]) < atmp) {x[i] = atmp; if (idx) idx[i] = j;} 1910 } 1911 } 1912 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 1913 PetscFunctionReturn(0); 1914 } 1915 1916 #undef __FUNCT__ 1917 #define __FUNCT__ "MatGetRowMin_SeqDense" 1918 PetscErrorCode MatGetRowMin_SeqDense(Mat A,Vec v,PetscInt idx[]) 1919 { 1920 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1921 PetscErrorCode ierr; 1922 PetscInt i,j,m = A->rmap->n,n = A->cmap->n,p; 1923 PetscScalar *x; 1924 MatScalar *aa = a->v; 1925 1926 PetscFunctionBegin; 1927 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 1928 1929 ierr = VecSet(v,0.0);CHKERRQ(ierr); 1930 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 1931 ierr = VecGetLocalSize(v,&p);CHKERRQ(ierr); 1932 if (p != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 1933 for (i=0; i<m; i++) { 1934 x[i] = aa[i]; if (idx) idx[i] = 0; 1935 for (j=1; j<n; j++) { 1936 if (PetscRealPart(x[i]) > PetscRealPart(aa[i+m*j])) {x[i] = aa[i + m*j]; if (idx) idx[i] = j;} 1937 } 1938 } 1939 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 1940 PetscFunctionReturn(0); 1941 } 1942 1943 #undef __FUNCT__ 1944 #define __FUNCT__ "MatGetColumnVector_SeqDense" 1945 PetscErrorCode MatGetColumnVector_SeqDense(Mat A,Vec v,PetscInt col) 1946 { 1947 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1948 PetscErrorCode ierr; 1949 PetscScalar *x; 1950 1951 PetscFunctionBegin; 1952 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 1953 1954 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 1955 ierr = PetscMemcpy(x,a->v+col*a->lda,A->rmap->n*sizeof(PetscScalar));CHKERRQ(ierr); 1956 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 1957 PetscFunctionReturn(0); 1958 } 1959 1960 1961 #undef __FUNCT__ 1962 #define __FUNCT__ "MatGetColumnNorms_SeqDense" 1963 PetscErrorCode MatGetColumnNorms_SeqDense(Mat A,NormType type,PetscReal *norms) 1964 { 1965 PetscErrorCode ierr; 1966 PetscInt i,j,m,n; 1967 PetscScalar *a; 1968 1969 PetscFunctionBegin; 1970 ierr = MatGetSize(A,&m,&n);CHKERRQ(ierr); 1971 ierr = PetscMemzero(norms,n*sizeof(PetscReal));CHKERRQ(ierr); 1972 ierr = MatDenseGetArray(A,&a);CHKERRQ(ierr); 1973 if (type == NORM_2) { 1974 for (i=0; i<n; i++) { 1975 for (j=0; j<m; j++) { 1976 norms[i] += PetscAbsScalar(a[j]*a[j]); 1977 } 1978 a += m; 1979 } 1980 } else if (type == NORM_1) { 1981 for (i=0; i<n; i++) { 1982 for (j=0; j<m; j++) { 1983 norms[i] += PetscAbsScalar(a[j]); 1984 } 1985 a += m; 1986 } 1987 } else if (type == NORM_INFINITY) { 1988 for (i=0; i<n; i++) { 1989 for (j=0; j<m; j++) { 1990 norms[i] = PetscMax(PetscAbsScalar(a[j]),norms[i]); 1991 } 1992 a += m; 1993 } 1994 } else SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Unknown NormType"); 1995 ierr = MatDenseRestoreArray(A,&a);CHKERRQ(ierr); 1996 if (type == NORM_2) { 1997 for (i=0; i<n; i++) norms[i] = PetscSqrtReal(norms[i]); 1998 } 1999 PetscFunctionReturn(0); 2000 } 2001 2002 #undef __FUNCT__ 2003 #define __FUNCT__ "MatSetRandom_SeqDense" 2004 static PetscErrorCode MatSetRandom_SeqDense(Mat x,PetscRandom rctx) 2005 { 2006 PetscErrorCode ierr; 2007 PetscScalar *a; 2008 PetscInt m,n,i; 2009 2010 PetscFunctionBegin; 2011 ierr = MatGetSize(x,&m,&n);CHKERRQ(ierr); 2012 ierr = MatDenseGetArray(x,&a);CHKERRQ(ierr); 2013 for (i=0; i<m*n; i++) { 2014 ierr = PetscRandomGetValue(rctx,a+i);CHKERRQ(ierr); 2015 } 2016 ierr = MatDenseRestoreArray(x,&a);CHKERRQ(ierr); 2017 PetscFunctionReturn(0); 2018 } 2019 2020 2021 /* -------------------------------------------------------------------*/ 2022 static struct _MatOps MatOps_Values = { MatSetValues_SeqDense, 2023 MatGetRow_SeqDense, 2024 MatRestoreRow_SeqDense, 2025 MatMult_SeqDense, 2026 /* 4*/ MatMultAdd_SeqDense, 2027 MatMultTranspose_SeqDense, 2028 MatMultTransposeAdd_SeqDense, 2029 0, 2030 0, 2031 0, 2032 /* 10*/ 0, 2033 MatLUFactor_SeqDense, 2034 MatCholeskyFactor_SeqDense, 2035 MatSOR_SeqDense, 2036 MatTranspose_SeqDense, 2037 /* 15*/ MatGetInfo_SeqDense, 2038 MatEqual_SeqDense, 2039 MatGetDiagonal_SeqDense, 2040 MatDiagonalScale_SeqDense, 2041 MatNorm_SeqDense, 2042 /* 20*/ MatAssemblyBegin_SeqDense, 2043 MatAssemblyEnd_SeqDense, 2044 MatSetOption_SeqDense, 2045 MatZeroEntries_SeqDense, 2046 /* 24*/ MatZeroRows_SeqDense, 2047 0, 2048 0, 2049 0, 2050 0, 2051 /* 29*/ MatSetUp_SeqDense, 2052 0, 2053 0, 2054 0, 2055 0, 2056 /* 34*/ MatDuplicate_SeqDense, 2057 0, 2058 0, 2059 0, 2060 0, 2061 /* 39*/ MatAXPY_SeqDense, 2062 MatGetSubMatrices_SeqDense, 2063 0, 2064 MatGetValues_SeqDense, 2065 MatCopy_SeqDense, 2066 /* 44*/ MatGetRowMax_SeqDense, 2067 MatScale_SeqDense, 2068 MatShift_Basic, 2069 0, 2070 0, 2071 /* 49*/ MatSetRandom_SeqDense, 2072 0, 2073 0, 2074 0, 2075 0, 2076 /* 54*/ 0, 2077 0, 2078 0, 2079 0, 2080 0, 2081 /* 59*/ 0, 2082 MatDestroy_SeqDense, 2083 MatView_SeqDense, 2084 0, 2085 0, 2086 /* 64*/ 0, 2087 0, 2088 0, 2089 0, 2090 0, 2091 /* 69*/ MatGetRowMaxAbs_SeqDense, 2092 0, 2093 0, 2094 0, 2095 0, 2096 /* 74*/ 0, 2097 0, 2098 0, 2099 0, 2100 0, 2101 /* 79*/ 0, 2102 0, 2103 0, 2104 0, 2105 /* 83*/ MatLoad_SeqDense, 2106 0, 2107 MatIsHermitian_SeqDense, 2108 0, 2109 0, 2110 0, 2111 /* 89*/ MatMatMult_SeqDense_SeqDense, 2112 MatMatMultSymbolic_SeqDense_SeqDense, 2113 MatMatMultNumeric_SeqDense_SeqDense, 2114 0, 2115 0, 2116 /* 94*/ 0, 2117 0, 2118 0, 2119 0, 2120 0, 2121 /* 99*/ 0, 2122 0, 2123 0, 2124 MatConjugate_SeqDense, 2125 0, 2126 /*104*/ 0, 2127 MatRealPart_SeqDense, 2128 MatImaginaryPart_SeqDense, 2129 0, 2130 0, 2131 /*109*/ MatMatSolve_SeqDense, 2132 0, 2133 MatGetRowMin_SeqDense, 2134 MatGetColumnVector_SeqDense, 2135 0, 2136 /*114*/ 0, 2137 0, 2138 0, 2139 0, 2140 0, 2141 /*119*/ 0, 2142 0, 2143 0, 2144 0, 2145 0, 2146 /*124*/ 0, 2147 MatGetColumnNorms_SeqDense, 2148 0, 2149 0, 2150 0, 2151 /*129*/ 0, 2152 MatTransposeMatMult_SeqDense_SeqDense, 2153 MatTransposeMatMultSymbolic_SeqDense_SeqDense, 2154 MatTransposeMatMultNumeric_SeqDense_SeqDense, 2155 0, 2156 /*134*/ 0, 2157 0, 2158 0, 2159 0, 2160 0, 2161 /*139*/ 0, 2162 0, 2163 0 2164 }; 2165 2166 #undef __FUNCT__ 2167 #define __FUNCT__ "MatCreateSeqDense" 2168 /*@C 2169 MatCreateSeqDense - Creates a sequential dense matrix that 2170 is stored in column major order (the usual Fortran 77 manner). Many 2171 of the matrix operations use the BLAS and LAPACK routines. 2172 2173 Collective on MPI_Comm 2174 2175 Input Parameters: 2176 + comm - MPI communicator, set to PETSC_COMM_SELF 2177 . m - number of rows 2178 . n - number of columns 2179 - data - optional location of matrix data in column major order. Set data=NULL for PETSc 2180 to control all matrix memory allocation. 2181 2182 Output Parameter: 2183 . A - the matrix 2184 2185 Notes: 2186 The data input variable is intended primarily for Fortran programmers 2187 who wish to allocate their own matrix memory space. Most users should 2188 set data=NULL. 2189 2190 Level: intermediate 2191 2192 .keywords: dense, matrix, LAPACK, BLAS 2193 2194 .seealso: MatCreate(), MatCreateDense(), MatSetValues() 2195 @*/ 2196 PetscErrorCode MatCreateSeqDense(MPI_Comm comm,PetscInt m,PetscInt n,PetscScalar *data,Mat *A) 2197 { 2198 PetscErrorCode ierr; 2199 2200 PetscFunctionBegin; 2201 ierr = MatCreate(comm,A);CHKERRQ(ierr); 2202 ierr = MatSetSizes(*A,m,n,m,n);CHKERRQ(ierr); 2203 ierr = MatSetType(*A,MATSEQDENSE);CHKERRQ(ierr); 2204 ierr = MatSeqDenseSetPreallocation(*A,data);CHKERRQ(ierr); 2205 PetscFunctionReturn(0); 2206 } 2207 2208 #undef __FUNCT__ 2209 #define __FUNCT__ "MatSeqDenseSetPreallocation" 2210 /*@C 2211 MatSeqDenseSetPreallocation - Sets the array used for storing the matrix elements 2212 2213 Collective on MPI_Comm 2214 2215 Input Parameters: 2216 + B - the matrix 2217 - data - the array (or NULL) 2218 2219 Notes: 2220 The data input variable is intended primarily for Fortran programmers 2221 who wish to allocate their own matrix memory space. Most users should 2222 need not call this routine. 2223 2224 Level: intermediate 2225 2226 .keywords: dense, matrix, LAPACK, BLAS 2227 2228 .seealso: MatCreate(), MatCreateDense(), MatSetValues(), MatSeqDenseSetLDA() 2229 2230 @*/ 2231 PetscErrorCode MatSeqDenseSetPreallocation(Mat B,PetscScalar data[]) 2232 { 2233 PetscErrorCode ierr; 2234 2235 PetscFunctionBegin; 2236 ierr = PetscTryMethod(B,"MatSeqDenseSetPreallocation_C",(Mat,PetscScalar[]),(B,data));CHKERRQ(ierr); 2237 PetscFunctionReturn(0); 2238 } 2239 2240 #undef __FUNCT__ 2241 #define __FUNCT__ "MatSeqDenseSetPreallocation_SeqDense" 2242 PetscErrorCode MatSeqDenseSetPreallocation_SeqDense(Mat B,PetscScalar *data) 2243 { 2244 Mat_SeqDense *b; 2245 PetscErrorCode ierr; 2246 2247 PetscFunctionBegin; 2248 B->preallocated = PETSC_TRUE; 2249 2250 ierr = PetscLayoutSetUp(B->rmap);CHKERRQ(ierr); 2251 ierr = PetscLayoutSetUp(B->cmap);CHKERRQ(ierr); 2252 2253 b = (Mat_SeqDense*)B->data; 2254 b->Mmax = B->rmap->n; 2255 b->Nmax = B->cmap->n; 2256 if (b->lda <= 0 || b->changelda) b->lda = B->rmap->n; 2257 2258 if (!data) { /* petsc-allocated storage */ 2259 if (!b->user_alloc) { ierr = PetscFree(b->v);CHKERRQ(ierr); } 2260 ierr = PetscCalloc1((size_t)b->lda*b->Nmax,&b->v);CHKERRQ(ierr); 2261 ierr = PetscLogObjectMemory((PetscObject)B,b->lda*b->Nmax*sizeof(PetscScalar));CHKERRQ(ierr); 2262 2263 b->user_alloc = PETSC_FALSE; 2264 } else { /* user-allocated storage */ 2265 if (!b->user_alloc) { ierr = PetscFree(b->v);CHKERRQ(ierr); } 2266 b->v = data; 2267 b->user_alloc = PETSC_TRUE; 2268 } 2269 B->assembled = PETSC_TRUE; 2270 PetscFunctionReturn(0); 2271 } 2272 2273 #if defined(PETSC_HAVE_ELEMENTAL) 2274 #undef __FUNCT__ 2275 #define __FUNCT__ "MatConvert_SeqDense_Elemental" 2276 PETSC_EXTERN PetscErrorCode MatConvert_SeqDense_Elemental(Mat A, MatType newtype,MatReuse reuse,Mat *newmat) 2277 { 2278 Mat mat_elemental; 2279 PetscErrorCode ierr; 2280 PetscScalar *array,*v_colwise; 2281 PetscInt M=A->rmap->N,N=A->cmap->N,i,j,k,*rows,*cols; 2282 2283 PetscFunctionBegin; 2284 ierr = PetscMalloc3(M*N,&v_colwise,M,&rows,N,&cols);CHKERRQ(ierr); 2285 ierr = MatDenseGetArray(A,&array);CHKERRQ(ierr); 2286 /* convert column-wise array into row-wise v_colwise, see MatSetValues_Elemental() */ 2287 k = 0; 2288 for (j=0; j<N; j++) { 2289 cols[j] = j; 2290 for (i=0; i<M; i++) { 2291 v_colwise[j*M+i] = array[k++]; 2292 } 2293 } 2294 for (i=0; i<M; i++) { 2295 rows[i] = i; 2296 } 2297 ierr = MatDenseRestoreArray(A,&array);CHKERRQ(ierr); 2298 2299 ierr = MatCreate(PetscObjectComm((PetscObject)A), &mat_elemental);CHKERRQ(ierr); 2300 ierr = MatSetSizes(mat_elemental,PETSC_DECIDE,PETSC_DECIDE,M,N);CHKERRQ(ierr); 2301 ierr = MatSetType(mat_elemental,MATELEMENTAL);CHKERRQ(ierr); 2302 ierr = MatSetUp(mat_elemental);CHKERRQ(ierr); 2303 2304 /* PETSc-Elemental interaface uses axpy for setting off-processor entries, only ADD_VALUES is allowed */ 2305 ierr = MatSetValues(mat_elemental,M,rows,N,cols,v_colwise,ADD_VALUES);CHKERRQ(ierr); 2306 ierr = MatAssemblyBegin(mat_elemental, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2307 ierr = MatAssemblyEnd(mat_elemental, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2308 ierr = PetscFree3(v_colwise,rows,cols);CHKERRQ(ierr); 2309 2310 if (reuse == MAT_REUSE_MATRIX) { 2311 ierr = MatHeaderReplace(A,mat_elemental);CHKERRQ(ierr); 2312 } else { 2313 *newmat = mat_elemental; 2314 } 2315 PetscFunctionReturn(0); 2316 } 2317 #endif 2318 2319 #undef __FUNCT__ 2320 #define __FUNCT__ "MatSeqDenseSetLDA" 2321 /*@C 2322 MatSeqDenseSetLDA - Declare the leading dimension of the user-provided array 2323 2324 Input parameter: 2325 + A - the matrix 2326 - lda - the leading dimension 2327 2328 Notes: 2329 This routine is to be used in conjunction with MatSeqDenseSetPreallocation(); 2330 it asserts that the preallocation has a leading dimension (the LDA parameter 2331 of Blas and Lapack fame) larger than M, the first dimension of the matrix. 2332 2333 Level: intermediate 2334 2335 .keywords: dense, matrix, LAPACK, BLAS 2336 2337 .seealso: MatCreate(), MatCreateSeqDense(), MatSeqDenseSetPreallocation(), MatSetMaximumSize() 2338 2339 @*/ 2340 PetscErrorCode MatSeqDenseSetLDA(Mat B,PetscInt lda) 2341 { 2342 Mat_SeqDense *b = (Mat_SeqDense*)B->data; 2343 2344 PetscFunctionBegin; 2345 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); 2346 b->lda = lda; 2347 b->changelda = PETSC_FALSE; 2348 b->Mmax = PetscMax(b->Mmax,lda); 2349 PetscFunctionReturn(0); 2350 } 2351 2352 /*MC 2353 MATSEQDENSE - MATSEQDENSE = "seqdense" - A matrix type to be used for sequential dense matrices. 2354 2355 Options Database Keys: 2356 . -mat_type seqdense - sets the matrix type to "seqdense" during a call to MatSetFromOptions() 2357 2358 Level: beginner 2359 2360 .seealso: MatCreateSeqDense() 2361 2362 M*/ 2363 2364 #undef __FUNCT__ 2365 #define __FUNCT__ "MatCreate_SeqDense" 2366 PETSC_EXTERN PetscErrorCode MatCreate_SeqDense(Mat B) 2367 { 2368 Mat_SeqDense *b; 2369 PetscErrorCode ierr; 2370 PetscMPIInt size; 2371 2372 PetscFunctionBegin; 2373 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)B),&size);CHKERRQ(ierr); 2374 if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Comm must be of size 1"); 2375 2376 ierr = PetscNewLog(B,&b);CHKERRQ(ierr); 2377 ierr = PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 2378 B->data = (void*)b; 2379 2380 b->pivots = 0; 2381 b->roworiented = PETSC_TRUE; 2382 b->v = 0; 2383 b->changelda = PETSC_FALSE; 2384 2385 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseGetArray_C",MatDenseGetArray_SeqDense);CHKERRQ(ierr); 2386 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreArray_C",MatDenseRestoreArray_SeqDense);CHKERRQ(ierr); 2387 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqdense_seqaij_C",MatConvert_SeqDense_SeqAIJ);CHKERRQ(ierr); 2388 #if defined(PETSC_HAVE_ELEMENTAL) 2389 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqdense_elemental_C",MatConvert_SeqDense_Elemental);CHKERRQ(ierr); 2390 #endif 2391 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqDenseSetPreallocation_C",MatSeqDenseSetPreallocation_SeqDense);CHKERRQ(ierr); 2392 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMult_seqaij_seqdense_C",MatMatMult_SeqAIJ_SeqDense);CHKERRQ(ierr); 2393 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMultSymbolic_seqaij_seqdense_C",MatMatMultSymbolic_SeqAIJ_SeqDense);CHKERRQ(ierr); 2394 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMultNumeric_seqaij_seqdense_C",MatMatMultNumeric_SeqAIJ_SeqDense);CHKERRQ(ierr); 2395 ierr = PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMult_seqaij_seqdense_C",MatTransposeMatMult_SeqAIJ_SeqDense);CHKERRQ(ierr); 2396 ierr = PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultSymbolic_seqaij_seqdense_C",MatTransposeMatMultSymbolic_SeqAIJ_SeqDense);CHKERRQ(ierr); 2397 ierr = PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultNumeric_seqaij_seqdense_C",MatTransposeMatMultNumeric_SeqAIJ_SeqDense);CHKERRQ(ierr); 2398 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQDENSE);CHKERRQ(ierr); 2399 PetscFunctionReturn(0); 2400 } 2401