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