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