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