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