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