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