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 (!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 a = (Mat_SeqDense*)newmat->data; 839 if (!a->user_alloc) { 840 ierr = MatSeqDenseSetPreallocation(newmat,NULL);CHKERRQ(ierr); 841 } 842 843 if (nz == MATRIX_BINARY_FORMAT_DENSE) { /* matrix in file is dense */ 844 a = (Mat_SeqDense*)newmat->data; 845 v = a->v; 846 /* Allocate some temp space to read in the values and then flip them 847 from row major to column major */ 848 ierr = PetscMalloc1((M*N > 0 ? M*N : 1),&w);CHKERRQ(ierr); 849 /* read in nonzero values */ 850 ierr = PetscBinaryRead(fd,w,M*N,PETSC_SCALAR);CHKERRQ(ierr); 851 /* now flip the values and store them in the matrix*/ 852 for (j=0; j<N; j++) { 853 for (i=0; i<M; i++) { 854 *v++ =w[i*N+j]; 855 } 856 } 857 ierr = PetscFree(w);CHKERRQ(ierr); 858 } else { 859 /* read row lengths */ 860 ierr = PetscMalloc1((M+1),&rowlengths);CHKERRQ(ierr); 861 ierr = PetscBinaryRead(fd,rowlengths,M,PETSC_INT);CHKERRQ(ierr); 862 863 a = (Mat_SeqDense*)newmat->data; 864 v = a->v; 865 866 /* read column indices and nonzeros */ 867 ierr = PetscMalloc1((nz+1),&scols);CHKERRQ(ierr); 868 cols = scols; 869 ierr = PetscBinaryRead(fd,cols,nz,PETSC_INT);CHKERRQ(ierr); 870 ierr = PetscMalloc1((nz+1),&svals);CHKERRQ(ierr); 871 vals = svals; 872 ierr = PetscBinaryRead(fd,vals,nz,PETSC_SCALAR);CHKERRQ(ierr); 873 874 /* insert into matrix */ 875 for (i=0; i<M; i++) { 876 for (j=0; j<rowlengths[i]; j++) v[i+M*scols[j]] = svals[j]; 877 svals += rowlengths[i]; scols += rowlengths[i]; 878 } 879 ierr = PetscFree(vals);CHKERRQ(ierr); 880 ierr = PetscFree(cols);CHKERRQ(ierr); 881 ierr = PetscFree(rowlengths);CHKERRQ(ierr); 882 } 883 ierr = MatAssemblyBegin(newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 884 ierr = MatAssemblyEnd(newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 885 PetscFunctionReturn(0); 886 } 887 888 #undef __FUNCT__ 889 #define __FUNCT__ "MatView_SeqDense_ASCII" 890 static PetscErrorCode MatView_SeqDense_ASCII(Mat A,PetscViewer viewer) 891 { 892 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 893 PetscErrorCode ierr; 894 PetscInt i,j; 895 const char *name; 896 PetscScalar *v; 897 PetscViewerFormat format; 898 #if defined(PETSC_USE_COMPLEX) 899 PetscBool allreal = PETSC_TRUE; 900 #endif 901 902 PetscFunctionBegin; 903 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 904 if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) { 905 PetscFunctionReturn(0); /* do nothing for now */ 906 } else if (format == PETSC_VIEWER_ASCII_COMMON) { 907 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);CHKERRQ(ierr); 908 for (i=0; i<A->rmap->n; i++) { 909 v = a->v + i; 910 ierr = PetscViewerASCIIPrintf(viewer,"row %D:",i);CHKERRQ(ierr); 911 for (j=0; j<A->cmap->n; j++) { 912 #if defined(PETSC_USE_COMPLEX) 913 if (PetscRealPart(*v) != 0.0 && PetscImaginaryPart(*v) != 0.0) { 914 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g + %g i) ",j,(double)PetscRealPart(*v),(double)PetscImaginaryPart(*v));CHKERRQ(ierr); 915 } else if (PetscRealPart(*v)) { 916 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g) ",j,(double)PetscRealPart(*v));CHKERRQ(ierr); 917 } 918 #else 919 if (*v) { 920 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g) ",j,(double)*v);CHKERRQ(ierr); 921 } 922 #endif 923 v += a->lda; 924 } 925 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 926 } 927 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);CHKERRQ(ierr); 928 } else { 929 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);CHKERRQ(ierr); 930 #if defined(PETSC_USE_COMPLEX) 931 /* determine if matrix has all real values */ 932 v = a->v; 933 for (i=0; i<A->rmap->n*A->cmap->n; i++) { 934 if (PetscImaginaryPart(v[i])) { allreal = PETSC_FALSE; break;} 935 } 936 #endif 937 if (format == PETSC_VIEWER_ASCII_MATLAB) { 938 ierr = PetscObjectGetName((PetscObject)A,&name);CHKERRQ(ierr); 939 ierr = PetscViewerASCIIPrintf(viewer,"%% Size = %D %D \n",A->rmap->n,A->cmap->n);CHKERRQ(ierr); 940 ierr = PetscViewerASCIIPrintf(viewer,"%s = zeros(%D,%D);\n",name,A->rmap->n,A->cmap->n);CHKERRQ(ierr); 941 ierr = PetscViewerASCIIPrintf(viewer,"%s = [\n",name);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,"MatConvert_seqdense_seqaij_C",NULL);CHKERRQ(ierr); 1169 #if defined(PETSC_HAVE_ELEMENTAL) 1170 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatConvert_seqdense_elemental_C",NULL);CHKERRQ(ierr); 1171 #endif 1172 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatGetFactor_petsc_C",NULL);CHKERRQ(ierr); 1173 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatSeqDenseSetPreallocation_C",NULL);CHKERRQ(ierr); 1174 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatMatMult_seqaij_seqdense_C",NULL);CHKERRQ(ierr); 1175 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatMatMultSymbolic_seqaij_seqdense_C",NULL);CHKERRQ(ierr); 1176 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatMatMultNumeric_seqaij_seqdense_C",NULL);CHKERRQ(ierr); 1177 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatTransposeMatMult_seqaij_seqdense_C",NULL);CHKERRQ(ierr); 1178 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatTransposeMatMultSymbolic_seqaij_seqdense_C",NULL);CHKERRQ(ierr); 1179 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatTransposeMatMultNumeric_seqaij_seqdense_C",NULL);CHKERRQ(ierr); 1180 PetscFunctionReturn(0); 1181 } 1182 1183 #undef __FUNCT__ 1184 #define __FUNCT__ "MatTranspose_SeqDense" 1185 PetscErrorCode MatTranspose_SeqDense(Mat A,MatReuse reuse,Mat *matout) 1186 { 1187 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1188 PetscErrorCode ierr; 1189 PetscInt k,j,m,n,M; 1190 PetscScalar *v,tmp; 1191 1192 PetscFunctionBegin; 1193 v = mat->v; m = A->rmap->n; M = mat->lda; n = A->cmap->n; 1194 if (reuse == MAT_REUSE_MATRIX && *matout == A) { /* in place transpose */ 1195 if (m != n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Can not transpose non-square matrix in place"); 1196 else { 1197 for (j=0; j<m; j++) { 1198 for (k=0; k<j; k++) { 1199 tmp = v[j + k*M]; 1200 v[j + k*M] = v[k + j*M]; 1201 v[k + j*M] = tmp; 1202 } 1203 } 1204 } 1205 } else { /* out-of-place transpose */ 1206 Mat tmat; 1207 Mat_SeqDense *tmatd; 1208 PetscScalar *v2; 1209 1210 if (reuse == MAT_INITIAL_MATRIX) { 1211 ierr = MatCreate(PetscObjectComm((PetscObject)A),&tmat);CHKERRQ(ierr); 1212 ierr = MatSetSizes(tmat,A->cmap->n,A->rmap->n,A->cmap->n,A->rmap->n);CHKERRQ(ierr); 1213 ierr = MatSetType(tmat,((PetscObject)A)->type_name);CHKERRQ(ierr); 1214 ierr = MatSeqDenseSetPreallocation(tmat,NULL);CHKERRQ(ierr); 1215 } else { 1216 tmat = *matout; 1217 } 1218 tmatd = (Mat_SeqDense*)tmat->data; 1219 v = mat->v; v2 = tmatd->v; 1220 for (j=0; j<n; j++) { 1221 for (k=0; k<m; k++) v2[j + k*n] = v[k + j*M]; 1222 } 1223 ierr = MatAssemblyBegin(tmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1224 ierr = MatAssemblyEnd(tmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1225 *matout = tmat; 1226 } 1227 PetscFunctionReturn(0); 1228 } 1229 1230 #undef __FUNCT__ 1231 #define __FUNCT__ "MatEqual_SeqDense" 1232 PetscErrorCode MatEqual_SeqDense(Mat A1,Mat A2,PetscBool *flg) 1233 { 1234 Mat_SeqDense *mat1 = (Mat_SeqDense*)A1->data; 1235 Mat_SeqDense *mat2 = (Mat_SeqDense*)A2->data; 1236 PetscInt i,j; 1237 PetscScalar *v1,*v2; 1238 1239 PetscFunctionBegin; 1240 if (A1->rmap->n != A2->rmap->n) {*flg = PETSC_FALSE; PetscFunctionReturn(0);} 1241 if (A1->cmap->n != A2->cmap->n) {*flg = PETSC_FALSE; PetscFunctionReturn(0);} 1242 for (i=0; i<A1->rmap->n; i++) { 1243 v1 = mat1->v+i; v2 = mat2->v+i; 1244 for (j=0; j<A1->cmap->n; j++) { 1245 if (*v1 != *v2) {*flg = PETSC_FALSE; PetscFunctionReturn(0);} 1246 v1 += mat1->lda; v2 += mat2->lda; 1247 } 1248 } 1249 *flg = PETSC_TRUE; 1250 PetscFunctionReturn(0); 1251 } 1252 1253 #undef __FUNCT__ 1254 #define __FUNCT__ "MatGetDiagonal_SeqDense" 1255 PetscErrorCode MatGetDiagonal_SeqDense(Mat A,Vec v) 1256 { 1257 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1258 PetscErrorCode ierr; 1259 PetscInt i,n,len; 1260 PetscScalar *x,zero = 0.0; 1261 1262 PetscFunctionBegin; 1263 ierr = VecSet(v,zero);CHKERRQ(ierr); 1264 ierr = VecGetSize(v,&n);CHKERRQ(ierr); 1265 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 1266 len = PetscMin(A->rmap->n,A->cmap->n); 1267 if (n != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming mat and vec"); 1268 for (i=0; i<len; i++) { 1269 x[i] = mat->v[i*mat->lda + i]; 1270 } 1271 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 1272 PetscFunctionReturn(0); 1273 } 1274 1275 #undef __FUNCT__ 1276 #define __FUNCT__ "MatDiagonalScale_SeqDense" 1277 PetscErrorCode MatDiagonalScale_SeqDense(Mat A,Vec ll,Vec rr) 1278 { 1279 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1280 PetscScalar *l,*r,x,*v; 1281 PetscErrorCode ierr; 1282 PetscInt i,j,m = A->rmap->n,n = A->cmap->n; 1283 1284 PetscFunctionBegin; 1285 if (ll) { 1286 ierr = VecGetSize(ll,&m);CHKERRQ(ierr); 1287 ierr = VecGetArray(ll,&l);CHKERRQ(ierr); 1288 if (m != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Left scaling vec wrong size"); 1289 for (i=0; i<m; i++) { 1290 x = l[i]; 1291 v = mat->v + i; 1292 for (j=0; j<n; j++) { (*v) *= x; v+= m;} 1293 } 1294 ierr = VecRestoreArray(ll,&l);CHKERRQ(ierr); 1295 ierr = PetscLogFlops(n*m);CHKERRQ(ierr); 1296 } 1297 if (rr) { 1298 ierr = VecGetSize(rr,&n);CHKERRQ(ierr); 1299 ierr = VecGetArray(rr,&r);CHKERRQ(ierr); 1300 if (n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Right scaling vec wrong size"); 1301 for (i=0; i<n; i++) { 1302 x = r[i]; 1303 v = mat->v + i*m; 1304 for (j=0; j<m; j++) (*v++) *= x; 1305 } 1306 ierr = VecRestoreArray(rr,&r);CHKERRQ(ierr); 1307 ierr = PetscLogFlops(n*m);CHKERRQ(ierr); 1308 } 1309 PetscFunctionReturn(0); 1310 } 1311 1312 #undef __FUNCT__ 1313 #define __FUNCT__ "MatNorm_SeqDense" 1314 PetscErrorCode MatNorm_SeqDense(Mat A,NormType type,PetscReal *nrm) 1315 { 1316 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1317 PetscScalar *v = mat->v; 1318 PetscReal sum = 0.0; 1319 PetscInt lda =mat->lda,m=A->rmap->n,i,j; 1320 PetscErrorCode ierr; 1321 1322 PetscFunctionBegin; 1323 if (type == NORM_FROBENIUS) { 1324 if (lda>m) { 1325 for (j=0; j<A->cmap->n; j++) { 1326 v = mat->v+j*lda; 1327 for (i=0; i<m; i++) { 1328 sum += PetscRealPart(PetscConj(*v)*(*v)); v++; 1329 } 1330 } 1331 } else { 1332 for (i=0; i<A->cmap->n*A->rmap->n; i++) { 1333 sum += PetscRealPart(PetscConj(*v)*(*v)); v++; 1334 } 1335 } 1336 *nrm = PetscSqrtReal(sum); 1337 ierr = PetscLogFlops(2.0*A->cmap->n*A->rmap->n);CHKERRQ(ierr); 1338 } else if (type == NORM_1) { 1339 *nrm = 0.0; 1340 for (j=0; j<A->cmap->n; j++) { 1341 v = mat->v + j*mat->lda; 1342 sum = 0.0; 1343 for (i=0; i<A->rmap->n; i++) { 1344 sum += PetscAbsScalar(*v); v++; 1345 } 1346 if (sum > *nrm) *nrm = sum; 1347 } 1348 ierr = PetscLogFlops(A->cmap->n*A->rmap->n);CHKERRQ(ierr); 1349 } else if (type == NORM_INFINITY) { 1350 *nrm = 0.0; 1351 for (j=0; j<A->rmap->n; j++) { 1352 v = mat->v + j; 1353 sum = 0.0; 1354 for (i=0; i<A->cmap->n; i++) { 1355 sum += PetscAbsScalar(*v); v += mat->lda; 1356 } 1357 if (sum > *nrm) *nrm = sum; 1358 } 1359 ierr = PetscLogFlops(A->cmap->n*A->rmap->n);CHKERRQ(ierr); 1360 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No two norm"); 1361 PetscFunctionReturn(0); 1362 } 1363 1364 #undef __FUNCT__ 1365 #define __FUNCT__ "MatSetOption_SeqDense" 1366 PetscErrorCode MatSetOption_SeqDense(Mat A,MatOption op,PetscBool flg) 1367 { 1368 Mat_SeqDense *aij = (Mat_SeqDense*)A->data; 1369 PetscErrorCode ierr; 1370 1371 PetscFunctionBegin; 1372 switch (op) { 1373 case MAT_ROW_ORIENTED: 1374 aij->roworiented = flg; 1375 break; 1376 case MAT_NEW_NONZERO_LOCATIONS: 1377 case MAT_NEW_NONZERO_LOCATION_ERR: 1378 case MAT_NEW_NONZERO_ALLOCATION_ERR: 1379 case MAT_NEW_DIAGONALS: 1380 case MAT_KEEP_NONZERO_PATTERN: 1381 case MAT_IGNORE_OFF_PROC_ENTRIES: 1382 case MAT_USE_HASH_TABLE: 1383 case MAT_IGNORE_LOWER_TRIANGULAR: 1384 ierr = PetscInfo1(A,"Option %s ignored\n",MatOptions[op]);CHKERRQ(ierr); 1385 break; 1386 case MAT_SPD: 1387 case MAT_SYMMETRIC: 1388 case MAT_STRUCTURALLY_SYMMETRIC: 1389 case MAT_HERMITIAN: 1390 case MAT_SYMMETRY_ETERNAL: 1391 /* These options are handled directly by MatSetOption() */ 1392 break; 1393 default: 1394 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"unknown option %s",MatOptions[op]); 1395 } 1396 PetscFunctionReturn(0); 1397 } 1398 1399 #undef __FUNCT__ 1400 #define __FUNCT__ "MatZeroEntries_SeqDense" 1401 PetscErrorCode MatZeroEntries_SeqDense(Mat A) 1402 { 1403 Mat_SeqDense *l = (Mat_SeqDense*)A->data; 1404 PetscErrorCode ierr; 1405 PetscInt lda=l->lda,m=A->rmap->n,j; 1406 1407 PetscFunctionBegin; 1408 if (lda>m) { 1409 for (j=0; j<A->cmap->n; j++) { 1410 ierr = PetscMemzero(l->v+j*lda,m*sizeof(PetscScalar));CHKERRQ(ierr); 1411 } 1412 } else { 1413 ierr = PetscMemzero(l->v,A->rmap->n*A->cmap->n*sizeof(PetscScalar));CHKERRQ(ierr); 1414 } 1415 PetscFunctionReturn(0); 1416 } 1417 1418 #undef __FUNCT__ 1419 #define __FUNCT__ "MatZeroRows_SeqDense" 1420 PetscErrorCode MatZeroRows_SeqDense(Mat A,PetscInt N,const PetscInt rows[],PetscScalar diag,Vec x,Vec b) 1421 { 1422 PetscErrorCode ierr; 1423 Mat_SeqDense *l = (Mat_SeqDense*)A->data; 1424 PetscInt m = l->lda, n = A->cmap->n, i,j; 1425 PetscScalar *slot,*bb; 1426 const PetscScalar *xx; 1427 1428 PetscFunctionBegin; 1429 #if defined(PETSC_USE_DEBUG) 1430 for (i=0; i<N; i++) { 1431 if (rows[i] < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative row requested to be zeroed"); 1432 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); 1433 } 1434 #endif 1435 1436 /* fix right hand side if needed */ 1437 if (x && b) { 1438 ierr = VecGetArrayRead(x,&xx);CHKERRQ(ierr); 1439 ierr = VecGetArray(b,&bb);CHKERRQ(ierr); 1440 for (i=0; i<N; i++) bb[rows[i]] = diag*xx[rows[i]]; 1441 ierr = VecRestoreArrayRead(x,&xx);CHKERRQ(ierr); 1442 ierr = VecRestoreArray(b,&bb);CHKERRQ(ierr); 1443 } 1444 1445 for (i=0; i<N; i++) { 1446 slot = l->v + rows[i]; 1447 for (j=0; j<n; j++) { *slot = 0.0; slot += m;} 1448 } 1449 if (diag != 0.0) { 1450 if (A->rmap->n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only coded for square matrices"); 1451 for (i=0; i<N; i++) { 1452 slot = l->v + (m+1)*rows[i]; 1453 *slot = diag; 1454 } 1455 } 1456 PetscFunctionReturn(0); 1457 } 1458 1459 #undef __FUNCT__ 1460 #define __FUNCT__ "MatDenseGetArray_SeqDense" 1461 PetscErrorCode MatDenseGetArray_SeqDense(Mat A,PetscScalar *array[]) 1462 { 1463 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1464 1465 PetscFunctionBegin; 1466 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"); 1467 *array = mat->v; 1468 PetscFunctionReturn(0); 1469 } 1470 1471 #undef __FUNCT__ 1472 #define __FUNCT__ "MatDenseRestoreArray_SeqDense" 1473 PetscErrorCode MatDenseRestoreArray_SeqDense(Mat A,PetscScalar *array[]) 1474 { 1475 PetscFunctionBegin; 1476 *array = 0; /* user cannot accidently use the array later */ 1477 PetscFunctionReturn(0); 1478 } 1479 1480 #undef __FUNCT__ 1481 #define __FUNCT__ "MatDenseGetArray" 1482 /*@C 1483 MatDenseGetArray - gives access to the array where the data for a SeqDense matrix is stored 1484 1485 Not Collective 1486 1487 Input Parameter: 1488 . mat - a MATSEQDENSE or MATMPIDENSE matrix 1489 1490 Output Parameter: 1491 . array - pointer to the data 1492 1493 Level: intermediate 1494 1495 .seealso: MatDenseRestoreArray() 1496 @*/ 1497 PetscErrorCode MatDenseGetArray(Mat A,PetscScalar **array) 1498 { 1499 PetscErrorCode ierr; 1500 1501 PetscFunctionBegin; 1502 ierr = PetscUseMethod(A,"MatDenseGetArray_C",(Mat,PetscScalar**),(A,array));CHKERRQ(ierr); 1503 PetscFunctionReturn(0); 1504 } 1505 1506 #undef __FUNCT__ 1507 #define __FUNCT__ "MatDenseRestoreArray" 1508 /*@C 1509 MatDenseRestoreArray - returns access to the array where the data for a dense matrix is stored obtained by MatDenseGetArray() 1510 1511 Not Collective 1512 1513 Input Parameters: 1514 . mat - a MATSEQDENSE or MATMPIDENSE matrix 1515 . array - pointer to the data 1516 1517 Level: intermediate 1518 1519 .seealso: MatDenseGetArray() 1520 @*/ 1521 PetscErrorCode MatDenseRestoreArray(Mat A,PetscScalar **array) 1522 { 1523 PetscErrorCode ierr; 1524 1525 PetscFunctionBegin; 1526 ierr = PetscUseMethod(A,"MatDenseRestoreArray_C",(Mat,PetscScalar**),(A,array));CHKERRQ(ierr); 1527 PetscFunctionReturn(0); 1528 } 1529 1530 #undef __FUNCT__ 1531 #define __FUNCT__ "MatGetSubMatrix_SeqDense" 1532 static PetscErrorCode MatGetSubMatrix_SeqDense(Mat A,IS isrow,IS iscol,PetscInt cs,MatReuse scall,Mat *B) 1533 { 1534 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1535 PetscErrorCode ierr; 1536 PetscInt i,j,nrows,ncols; 1537 const PetscInt *irow,*icol; 1538 PetscScalar *av,*bv,*v = mat->v; 1539 Mat newmat; 1540 1541 PetscFunctionBegin; 1542 ierr = ISGetIndices(isrow,&irow);CHKERRQ(ierr); 1543 ierr = ISGetIndices(iscol,&icol);CHKERRQ(ierr); 1544 ierr = ISGetLocalSize(isrow,&nrows);CHKERRQ(ierr); 1545 ierr = ISGetLocalSize(iscol,&ncols);CHKERRQ(ierr); 1546 1547 /* Check submatrixcall */ 1548 if (scall == MAT_REUSE_MATRIX) { 1549 PetscInt n_cols,n_rows; 1550 ierr = MatGetSize(*B,&n_rows,&n_cols);CHKERRQ(ierr); 1551 if (n_rows != nrows || n_cols != ncols) { 1552 /* resize the result matrix to match number of requested rows/columns */ 1553 ierr = MatSetSizes(*B,nrows,ncols,nrows,ncols);CHKERRQ(ierr); 1554 } 1555 newmat = *B; 1556 } else { 1557 /* Create and fill new matrix */ 1558 ierr = MatCreate(PetscObjectComm((PetscObject)A),&newmat);CHKERRQ(ierr); 1559 ierr = MatSetSizes(newmat,nrows,ncols,nrows,ncols);CHKERRQ(ierr); 1560 ierr = MatSetType(newmat,((PetscObject)A)->type_name);CHKERRQ(ierr); 1561 ierr = MatSeqDenseSetPreallocation(newmat,NULL);CHKERRQ(ierr); 1562 } 1563 1564 /* Now extract the data pointers and do the copy,column at a time */ 1565 bv = ((Mat_SeqDense*)newmat->data)->v; 1566 1567 for (i=0; i<ncols; i++) { 1568 av = v + mat->lda*icol[i]; 1569 for (j=0; j<nrows; j++) *bv++ = av[irow[j]]; 1570 } 1571 1572 /* Assemble the matrices so that the correct flags are set */ 1573 ierr = MatAssemblyBegin(newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1574 ierr = MatAssemblyEnd(newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1575 1576 /* Free work space */ 1577 ierr = ISRestoreIndices(isrow,&irow);CHKERRQ(ierr); 1578 ierr = ISRestoreIndices(iscol,&icol);CHKERRQ(ierr); 1579 *B = newmat; 1580 PetscFunctionReturn(0); 1581 } 1582 1583 #undef __FUNCT__ 1584 #define __FUNCT__ "MatGetSubMatrices_SeqDense" 1585 PetscErrorCode MatGetSubMatrices_SeqDense(Mat A,PetscInt n,const IS irow[],const IS icol[],MatReuse scall,Mat *B[]) 1586 { 1587 PetscErrorCode ierr; 1588 PetscInt i; 1589 1590 PetscFunctionBegin; 1591 if (scall == MAT_INITIAL_MATRIX) { 1592 ierr = PetscMalloc1((n+1),B);CHKERRQ(ierr); 1593 } 1594 1595 for (i=0; i<n; i++) { 1596 ierr = MatGetSubMatrix_SeqDense(A,irow[i],icol[i],PETSC_DECIDE,scall,&(*B)[i]);CHKERRQ(ierr); 1597 } 1598 PetscFunctionReturn(0); 1599 } 1600 1601 #undef __FUNCT__ 1602 #define __FUNCT__ "MatAssemblyBegin_SeqDense" 1603 PetscErrorCode MatAssemblyBegin_SeqDense(Mat mat,MatAssemblyType mode) 1604 { 1605 PetscFunctionBegin; 1606 PetscFunctionReturn(0); 1607 } 1608 1609 #undef __FUNCT__ 1610 #define __FUNCT__ "MatAssemblyEnd_SeqDense" 1611 PetscErrorCode MatAssemblyEnd_SeqDense(Mat mat,MatAssemblyType mode) 1612 { 1613 PetscFunctionBegin; 1614 PetscFunctionReturn(0); 1615 } 1616 1617 #undef __FUNCT__ 1618 #define __FUNCT__ "MatCopy_SeqDense" 1619 PetscErrorCode MatCopy_SeqDense(Mat A,Mat B,MatStructure str) 1620 { 1621 Mat_SeqDense *a = (Mat_SeqDense*)A->data,*b = (Mat_SeqDense*)B->data; 1622 PetscErrorCode ierr; 1623 PetscInt lda1=a->lda,lda2=b->lda, m=A->rmap->n,n=A->cmap->n, j; 1624 1625 PetscFunctionBegin; 1626 /* If the two matrices don't have the same copy implementation, they aren't compatible for fast copy. */ 1627 if (A->ops->copy != B->ops->copy) { 1628 ierr = MatCopy_Basic(A,B,str);CHKERRQ(ierr); 1629 PetscFunctionReturn(0); 1630 } 1631 if (m != B->rmap->n || n != B->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"size(B) != size(A)"); 1632 if (lda1>m || lda2>m) { 1633 for (j=0; j<n; j++) { 1634 ierr = PetscMemcpy(b->v+j*lda2,a->v+j*lda1,m*sizeof(PetscScalar));CHKERRQ(ierr); 1635 } 1636 } else { 1637 ierr = PetscMemcpy(b->v,a->v,A->rmap->n*A->cmap->n*sizeof(PetscScalar));CHKERRQ(ierr); 1638 } 1639 PetscFunctionReturn(0); 1640 } 1641 1642 #undef __FUNCT__ 1643 #define __FUNCT__ "MatSetUp_SeqDense" 1644 PetscErrorCode MatSetUp_SeqDense(Mat A) 1645 { 1646 PetscErrorCode ierr; 1647 1648 PetscFunctionBegin; 1649 ierr = MatSeqDenseSetPreallocation(A,0);CHKERRQ(ierr); 1650 PetscFunctionReturn(0); 1651 } 1652 1653 #undef __FUNCT__ 1654 #define __FUNCT__ "MatConjugate_SeqDense" 1655 static PetscErrorCode MatConjugate_SeqDense(Mat A) 1656 { 1657 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1658 PetscInt i,nz = A->rmap->n*A->cmap->n; 1659 PetscScalar *aa = a->v; 1660 1661 PetscFunctionBegin; 1662 for (i=0; i<nz; i++) aa[i] = PetscConj(aa[i]); 1663 PetscFunctionReturn(0); 1664 } 1665 1666 #undef __FUNCT__ 1667 #define __FUNCT__ "MatRealPart_SeqDense" 1668 static PetscErrorCode MatRealPart_SeqDense(Mat A) 1669 { 1670 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1671 PetscInt i,nz = A->rmap->n*A->cmap->n; 1672 PetscScalar *aa = a->v; 1673 1674 PetscFunctionBegin; 1675 for (i=0; i<nz; i++) aa[i] = PetscRealPart(aa[i]); 1676 PetscFunctionReturn(0); 1677 } 1678 1679 #undef __FUNCT__ 1680 #define __FUNCT__ "MatImaginaryPart_SeqDense" 1681 static PetscErrorCode MatImaginaryPart_SeqDense(Mat A) 1682 { 1683 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1684 PetscInt i,nz = A->rmap->n*A->cmap->n; 1685 PetscScalar *aa = a->v; 1686 1687 PetscFunctionBegin; 1688 for (i=0; i<nz; i++) aa[i] = PetscImaginaryPart(aa[i]); 1689 PetscFunctionReturn(0); 1690 } 1691 1692 /* ----------------------------------------------------------------*/ 1693 #undef __FUNCT__ 1694 #define __FUNCT__ "MatMatMult_SeqDense_SeqDense" 1695 PetscErrorCode MatMatMult_SeqDense_SeqDense(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C) 1696 { 1697 PetscErrorCode ierr; 1698 1699 PetscFunctionBegin; 1700 if (scall == MAT_INITIAL_MATRIX) { 1701 ierr = PetscLogEventBegin(MAT_MatMultSymbolic,A,B,0,0);CHKERRQ(ierr); 1702 ierr = MatMatMultSymbolic_SeqDense_SeqDense(A,B,fill,C);CHKERRQ(ierr); 1703 ierr = PetscLogEventEnd(MAT_MatMultSymbolic,A,B,0,0);CHKERRQ(ierr); 1704 } 1705 ierr = PetscLogEventBegin(MAT_MatMultNumeric,A,B,0,0);CHKERRQ(ierr); 1706 ierr = MatMatMultNumeric_SeqDense_SeqDense(A,B,*C);CHKERRQ(ierr); 1707 ierr = PetscLogEventEnd(MAT_MatMultNumeric,A,B,0,0);CHKERRQ(ierr); 1708 PetscFunctionReturn(0); 1709 } 1710 1711 #undef __FUNCT__ 1712 #define __FUNCT__ "MatMatMultSymbolic_SeqDense_SeqDense" 1713 PetscErrorCode MatMatMultSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat *C) 1714 { 1715 PetscErrorCode ierr; 1716 PetscInt m=A->rmap->n,n=B->cmap->n; 1717 Mat Cmat; 1718 1719 PetscFunctionBegin; 1720 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); 1721 ierr = MatCreate(PETSC_COMM_SELF,&Cmat);CHKERRQ(ierr); 1722 ierr = MatSetSizes(Cmat,m,n,m,n);CHKERRQ(ierr); 1723 ierr = MatSetType(Cmat,MATSEQDENSE);CHKERRQ(ierr); 1724 ierr = MatSeqDenseSetPreallocation(Cmat,NULL);CHKERRQ(ierr); 1725 1726 *C = Cmat; 1727 PetscFunctionReturn(0); 1728 } 1729 1730 #undef __FUNCT__ 1731 #define __FUNCT__ "MatMatMultNumeric_SeqDense_SeqDense" 1732 PetscErrorCode MatMatMultNumeric_SeqDense_SeqDense(Mat A,Mat B,Mat C) 1733 { 1734 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1735 Mat_SeqDense *b = (Mat_SeqDense*)B->data; 1736 Mat_SeqDense *c = (Mat_SeqDense*)C->data; 1737 PetscBLASInt m,n,k; 1738 PetscScalar _DOne=1.0,_DZero=0.0; 1739 PetscErrorCode ierr; 1740 1741 PetscFunctionBegin; 1742 ierr = PetscBLASIntCast(A->rmap->n,&m);CHKERRQ(ierr); 1743 ierr = PetscBLASIntCast(B->cmap->n,&n);CHKERRQ(ierr); 1744 ierr = PetscBLASIntCast(A->cmap->n,&k);CHKERRQ(ierr); 1745 PetscStackCallBLAS("BLASgemv",BLASgemm_("N","N",&m,&n,&k,&_DOne,a->v,&a->lda,b->v,&b->lda,&_DZero,c->v,&c->lda)); 1746 PetscFunctionReturn(0); 1747 } 1748 1749 #undef __FUNCT__ 1750 #define __FUNCT__ "MatTransposeMatMult_SeqDense_SeqDense" 1751 PetscErrorCode MatTransposeMatMult_SeqDense_SeqDense(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C) 1752 { 1753 PetscErrorCode ierr; 1754 1755 PetscFunctionBegin; 1756 if (scall == MAT_INITIAL_MATRIX) { 1757 ierr = PetscLogEventBegin(MAT_TransposeMatMultSymbolic,A,B,0,0);CHKERRQ(ierr); 1758 ierr = MatTransposeMatMultSymbolic_SeqDense_SeqDense(A,B,fill,C);CHKERRQ(ierr); 1759 ierr = PetscLogEventEnd(MAT_TransposeMatMultSymbolic,A,B,0,0);CHKERRQ(ierr); 1760 } 1761 ierr = PetscLogEventBegin(MAT_TransposeMatMultNumeric,A,B,0,0);CHKERRQ(ierr); 1762 ierr = MatTransposeMatMultNumeric_SeqDense_SeqDense(A,B,*C);CHKERRQ(ierr); 1763 ierr = PetscLogEventEnd(MAT_TransposeMatMultNumeric,A,B,0,0);CHKERRQ(ierr); 1764 PetscFunctionReturn(0); 1765 } 1766 1767 #undef __FUNCT__ 1768 #define __FUNCT__ "MatTransposeMatMultSymbolic_SeqDense_SeqDense" 1769 PetscErrorCode MatTransposeMatMultSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat *C) 1770 { 1771 PetscErrorCode ierr; 1772 PetscInt m=A->cmap->n,n=B->cmap->n; 1773 Mat Cmat; 1774 1775 PetscFunctionBegin; 1776 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); 1777 ierr = MatCreate(PETSC_COMM_SELF,&Cmat);CHKERRQ(ierr); 1778 ierr = MatSetSizes(Cmat,m,n,m,n);CHKERRQ(ierr); 1779 ierr = MatSetType(Cmat,MATSEQDENSE);CHKERRQ(ierr); 1780 ierr = MatSeqDenseSetPreallocation(Cmat,NULL);CHKERRQ(ierr); 1781 1782 Cmat->assembled = PETSC_TRUE; 1783 1784 *C = Cmat; 1785 PetscFunctionReturn(0); 1786 } 1787 1788 #undef __FUNCT__ 1789 #define __FUNCT__ "MatTransposeMatMultNumeric_SeqDense_SeqDense" 1790 PetscErrorCode MatTransposeMatMultNumeric_SeqDense_SeqDense(Mat A,Mat B,Mat C) 1791 { 1792 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1793 Mat_SeqDense *b = (Mat_SeqDense*)B->data; 1794 Mat_SeqDense *c = (Mat_SeqDense*)C->data; 1795 PetscBLASInt m,n,k; 1796 PetscScalar _DOne=1.0,_DZero=0.0; 1797 PetscErrorCode ierr; 1798 1799 PetscFunctionBegin; 1800 ierr = PetscBLASIntCast(A->cmap->n,&m);CHKERRQ(ierr); 1801 ierr = PetscBLASIntCast(B->cmap->n,&n);CHKERRQ(ierr); 1802 ierr = PetscBLASIntCast(A->rmap->n,&k);CHKERRQ(ierr); 1803 /* 1804 Note the m and n arguments below are the number rows and columns of A', not A! 1805 */ 1806 PetscStackCallBLAS("BLASgemv",BLASgemm_("T","N",&m,&n,&k,&_DOne,a->v,&a->lda,b->v,&b->lda,&_DZero,c->v,&c->lda)); 1807 PetscFunctionReturn(0); 1808 } 1809 1810 #undef __FUNCT__ 1811 #define __FUNCT__ "MatGetRowMax_SeqDense" 1812 PetscErrorCode MatGetRowMax_SeqDense(Mat A,Vec v,PetscInt idx[]) 1813 { 1814 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1815 PetscErrorCode ierr; 1816 PetscInt i,j,m = A->rmap->n,n = A->cmap->n,p; 1817 PetscScalar *x; 1818 MatScalar *aa = a->v; 1819 1820 PetscFunctionBegin; 1821 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 1822 1823 ierr = VecSet(v,0.0);CHKERRQ(ierr); 1824 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 1825 ierr = VecGetLocalSize(v,&p);CHKERRQ(ierr); 1826 if (p != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 1827 for (i=0; i<m; i++) { 1828 x[i] = aa[i]; if (idx) idx[i] = 0; 1829 for (j=1; j<n; j++) { 1830 if (PetscRealPart(x[i]) < PetscRealPart(aa[i+m*j])) {x[i] = aa[i + m*j]; if (idx) idx[i] = j;} 1831 } 1832 } 1833 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 1834 PetscFunctionReturn(0); 1835 } 1836 1837 #undef __FUNCT__ 1838 #define __FUNCT__ "MatGetRowMaxAbs_SeqDense" 1839 PetscErrorCode MatGetRowMaxAbs_SeqDense(Mat A,Vec v,PetscInt idx[]) 1840 { 1841 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1842 PetscErrorCode ierr; 1843 PetscInt i,j,m = A->rmap->n,n = A->cmap->n,p; 1844 PetscScalar *x; 1845 PetscReal atmp; 1846 MatScalar *aa = a->v; 1847 1848 PetscFunctionBegin; 1849 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 1850 1851 ierr = VecSet(v,0.0);CHKERRQ(ierr); 1852 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 1853 ierr = VecGetLocalSize(v,&p);CHKERRQ(ierr); 1854 if (p != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 1855 for (i=0; i<m; i++) { 1856 x[i] = PetscAbsScalar(aa[i]); 1857 for (j=1; j<n; j++) { 1858 atmp = PetscAbsScalar(aa[i+m*j]); 1859 if (PetscAbsScalar(x[i]) < atmp) {x[i] = atmp; if (idx) idx[i] = j;} 1860 } 1861 } 1862 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 1863 PetscFunctionReturn(0); 1864 } 1865 1866 #undef __FUNCT__ 1867 #define __FUNCT__ "MatGetRowMin_SeqDense" 1868 PetscErrorCode MatGetRowMin_SeqDense(Mat A,Vec v,PetscInt idx[]) 1869 { 1870 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1871 PetscErrorCode ierr; 1872 PetscInt i,j,m = A->rmap->n,n = A->cmap->n,p; 1873 PetscScalar *x; 1874 MatScalar *aa = a->v; 1875 1876 PetscFunctionBegin; 1877 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 1878 1879 ierr = VecSet(v,0.0);CHKERRQ(ierr); 1880 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 1881 ierr = VecGetLocalSize(v,&p);CHKERRQ(ierr); 1882 if (p != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 1883 for (i=0; i<m; i++) { 1884 x[i] = aa[i]; if (idx) idx[i] = 0; 1885 for (j=1; j<n; j++) { 1886 if (PetscRealPart(x[i]) > PetscRealPart(aa[i+m*j])) {x[i] = aa[i + m*j]; if (idx) idx[i] = j;} 1887 } 1888 } 1889 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 1890 PetscFunctionReturn(0); 1891 } 1892 1893 #undef __FUNCT__ 1894 #define __FUNCT__ "MatGetColumnVector_SeqDense" 1895 PetscErrorCode MatGetColumnVector_SeqDense(Mat A,Vec v,PetscInt col) 1896 { 1897 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1898 PetscErrorCode ierr; 1899 PetscScalar *x; 1900 1901 PetscFunctionBegin; 1902 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 1903 1904 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 1905 ierr = PetscMemcpy(x,a->v+col*a->lda,A->rmap->n*sizeof(PetscScalar));CHKERRQ(ierr); 1906 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 1907 PetscFunctionReturn(0); 1908 } 1909 1910 1911 #undef __FUNCT__ 1912 #define __FUNCT__ "MatGetColumnNorms_SeqDense" 1913 PetscErrorCode MatGetColumnNorms_SeqDense(Mat A,NormType type,PetscReal *norms) 1914 { 1915 PetscErrorCode ierr; 1916 PetscInt i,j,m,n; 1917 PetscScalar *a; 1918 1919 PetscFunctionBegin; 1920 ierr = MatGetSize(A,&m,&n);CHKERRQ(ierr); 1921 ierr = PetscMemzero(norms,n*sizeof(PetscReal));CHKERRQ(ierr); 1922 ierr = MatDenseGetArray(A,&a);CHKERRQ(ierr); 1923 if (type == NORM_2) { 1924 for (i=0; i<n; i++) { 1925 for (j=0; j<m; j++) { 1926 norms[i] += PetscAbsScalar(a[j]*a[j]); 1927 } 1928 a += m; 1929 } 1930 } else if (type == NORM_1) { 1931 for (i=0; i<n; i++) { 1932 for (j=0; j<m; j++) { 1933 norms[i] += PetscAbsScalar(a[j]); 1934 } 1935 a += m; 1936 } 1937 } else if (type == NORM_INFINITY) { 1938 for (i=0; i<n; i++) { 1939 for (j=0; j<m; j++) { 1940 norms[i] = PetscMax(PetscAbsScalar(a[j]),norms[i]); 1941 } 1942 a += m; 1943 } 1944 } else SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Unknown NormType"); 1945 ierr = MatDenseRestoreArray(A,&a);CHKERRQ(ierr); 1946 if (type == NORM_2) { 1947 for (i=0; i<n; i++) norms[i] = PetscSqrtReal(norms[i]); 1948 } 1949 PetscFunctionReturn(0); 1950 } 1951 1952 #undef __FUNCT__ 1953 #define __FUNCT__ "MatSetRandom_SeqDense" 1954 static PetscErrorCode MatSetRandom_SeqDense(Mat x,PetscRandom rctx) 1955 { 1956 PetscErrorCode ierr; 1957 PetscScalar *a; 1958 PetscInt m,n,i; 1959 1960 PetscFunctionBegin; 1961 ierr = MatGetSize(x,&m,&n);CHKERRQ(ierr); 1962 ierr = MatDenseGetArray(x,&a);CHKERRQ(ierr); 1963 for (i=0; i<m*n; i++) { 1964 ierr = PetscRandomGetValue(rctx,a+i);CHKERRQ(ierr); 1965 } 1966 ierr = MatDenseRestoreArray(x,&a);CHKERRQ(ierr); 1967 PetscFunctionReturn(0); 1968 } 1969 1970 1971 /* -------------------------------------------------------------------*/ 1972 static struct _MatOps MatOps_Values = { MatSetValues_SeqDense, 1973 MatGetRow_SeqDense, 1974 MatRestoreRow_SeqDense, 1975 MatMult_SeqDense, 1976 /* 4*/ MatMultAdd_SeqDense, 1977 MatMultTranspose_SeqDense, 1978 MatMultTransposeAdd_SeqDense, 1979 0, 1980 0, 1981 0, 1982 /* 10*/ 0, 1983 MatLUFactor_SeqDense, 1984 MatCholeskyFactor_SeqDense, 1985 MatSOR_SeqDense, 1986 MatTranspose_SeqDense, 1987 /* 15*/ MatGetInfo_SeqDense, 1988 MatEqual_SeqDense, 1989 MatGetDiagonal_SeqDense, 1990 MatDiagonalScale_SeqDense, 1991 MatNorm_SeqDense, 1992 /* 20*/ MatAssemblyBegin_SeqDense, 1993 MatAssemblyEnd_SeqDense, 1994 MatSetOption_SeqDense, 1995 MatZeroEntries_SeqDense, 1996 /* 24*/ MatZeroRows_SeqDense, 1997 0, 1998 0, 1999 0, 2000 0, 2001 /* 29*/ MatSetUp_SeqDense, 2002 0, 2003 0, 2004 0, 2005 0, 2006 /* 34*/ MatDuplicate_SeqDense, 2007 0, 2008 0, 2009 0, 2010 0, 2011 /* 39*/ MatAXPY_SeqDense, 2012 MatGetSubMatrices_SeqDense, 2013 0, 2014 MatGetValues_SeqDense, 2015 MatCopy_SeqDense, 2016 /* 44*/ MatGetRowMax_SeqDense, 2017 MatScale_SeqDense, 2018 0, 2019 0, 2020 0, 2021 /* 49*/ MatSetRandom_SeqDense, 2022 0, 2023 0, 2024 0, 2025 0, 2026 /* 54*/ 0, 2027 0, 2028 0, 2029 0, 2030 0, 2031 /* 59*/ 0, 2032 MatDestroy_SeqDense, 2033 MatView_SeqDense, 2034 0, 2035 0, 2036 /* 64*/ 0, 2037 0, 2038 0, 2039 0, 2040 0, 2041 /* 69*/ MatGetRowMaxAbs_SeqDense, 2042 0, 2043 0, 2044 0, 2045 0, 2046 /* 74*/ 0, 2047 0, 2048 0, 2049 0, 2050 0, 2051 /* 79*/ 0, 2052 0, 2053 0, 2054 0, 2055 /* 83*/ MatLoad_SeqDense, 2056 0, 2057 MatIsHermitian_SeqDense, 2058 0, 2059 0, 2060 0, 2061 /* 89*/ MatMatMult_SeqDense_SeqDense, 2062 MatMatMultSymbolic_SeqDense_SeqDense, 2063 MatMatMultNumeric_SeqDense_SeqDense, 2064 0, 2065 0, 2066 /* 94*/ 0, 2067 0, 2068 0, 2069 0, 2070 0, 2071 /* 99*/ 0, 2072 0, 2073 0, 2074 MatConjugate_SeqDense, 2075 0, 2076 /*104*/ 0, 2077 MatRealPart_SeqDense, 2078 MatImaginaryPart_SeqDense, 2079 0, 2080 0, 2081 /*109*/ MatMatSolve_SeqDense, 2082 0, 2083 MatGetRowMin_SeqDense, 2084 MatGetColumnVector_SeqDense, 2085 0, 2086 /*114*/ 0, 2087 0, 2088 0, 2089 0, 2090 0, 2091 /*119*/ 0, 2092 0, 2093 0, 2094 0, 2095 0, 2096 /*124*/ 0, 2097 MatGetColumnNorms_SeqDense, 2098 0, 2099 0, 2100 0, 2101 /*129*/ 0, 2102 MatTransposeMatMult_SeqDense_SeqDense, 2103 MatTransposeMatMultSymbolic_SeqDense_SeqDense, 2104 MatTransposeMatMultNumeric_SeqDense_SeqDense, 2105 0, 2106 /*134*/ 0, 2107 0, 2108 0, 2109 0, 2110 0, 2111 /*139*/ 0, 2112 0, 2113 0 2114 }; 2115 2116 #undef __FUNCT__ 2117 #define __FUNCT__ "MatCreateSeqDense" 2118 /*@C 2119 MatCreateSeqDense - Creates a sequential dense matrix that 2120 is stored in column major order (the usual Fortran 77 manner). Many 2121 of the matrix operations use the BLAS and LAPACK routines. 2122 2123 Collective on MPI_Comm 2124 2125 Input Parameters: 2126 + comm - MPI communicator, set to PETSC_COMM_SELF 2127 . m - number of rows 2128 . n - number of columns 2129 - data - optional location of matrix data in column major order. Set data=NULL for PETSc 2130 to control all matrix memory allocation. 2131 2132 Output Parameter: 2133 . A - the matrix 2134 2135 Notes: 2136 The data input variable is intended primarily for Fortran programmers 2137 who wish to allocate their own matrix memory space. Most users should 2138 set data=NULL. 2139 2140 Level: intermediate 2141 2142 .keywords: dense, matrix, LAPACK, BLAS 2143 2144 .seealso: MatCreate(), MatCreateDense(), MatSetValues() 2145 @*/ 2146 PetscErrorCode MatCreateSeqDense(MPI_Comm comm,PetscInt m,PetscInt n,PetscScalar *data,Mat *A) 2147 { 2148 PetscErrorCode ierr; 2149 2150 PetscFunctionBegin; 2151 ierr = MatCreate(comm,A);CHKERRQ(ierr); 2152 ierr = MatSetSizes(*A,m,n,m,n);CHKERRQ(ierr); 2153 ierr = MatSetType(*A,MATSEQDENSE);CHKERRQ(ierr); 2154 ierr = MatSeqDenseSetPreallocation(*A,data);CHKERRQ(ierr); 2155 PetscFunctionReturn(0); 2156 } 2157 2158 #undef __FUNCT__ 2159 #define __FUNCT__ "MatSeqDenseSetPreallocation" 2160 /*@C 2161 MatSeqDenseSetPreallocation - Sets the array used for storing the matrix elements 2162 2163 Collective on MPI_Comm 2164 2165 Input Parameters: 2166 + B - the matrix 2167 - data - the array (or NULL) 2168 2169 Notes: 2170 The data input variable is intended primarily for Fortran programmers 2171 who wish to allocate their own matrix memory space. Most users should 2172 need not call this routine. 2173 2174 Level: intermediate 2175 2176 .keywords: dense, matrix, LAPACK, BLAS 2177 2178 .seealso: MatCreate(), MatCreateDense(), MatSetValues(), MatSeqDenseSetLDA() 2179 2180 @*/ 2181 PetscErrorCode MatSeqDenseSetPreallocation(Mat B,PetscScalar data[]) 2182 { 2183 PetscErrorCode ierr; 2184 2185 PetscFunctionBegin; 2186 ierr = PetscTryMethod(B,"MatSeqDenseSetPreallocation_C",(Mat,PetscScalar[]),(B,data));CHKERRQ(ierr); 2187 PetscFunctionReturn(0); 2188 } 2189 2190 #undef __FUNCT__ 2191 #define __FUNCT__ "MatSeqDenseSetPreallocation_SeqDense" 2192 PetscErrorCode MatSeqDenseSetPreallocation_SeqDense(Mat B,PetscScalar *data) 2193 { 2194 Mat_SeqDense *b; 2195 PetscErrorCode ierr; 2196 2197 PetscFunctionBegin; 2198 B->preallocated = PETSC_TRUE; 2199 2200 ierr = PetscLayoutSetUp(B->rmap);CHKERRQ(ierr); 2201 ierr = PetscLayoutSetUp(B->cmap);CHKERRQ(ierr); 2202 2203 b = (Mat_SeqDense*)B->data; 2204 b->Mmax = B->rmap->n; 2205 b->Nmax = B->cmap->n; 2206 if (b->lda <= 0 || b->changelda) b->lda = B->rmap->n; 2207 2208 if (!data) { /* petsc-allocated storage */ 2209 if (!b->user_alloc) { ierr = PetscFree(b->v);CHKERRQ(ierr); } 2210 ierr = PetscCalloc1((size_t)b->lda*b->Nmax,&b->v);CHKERRQ(ierr); 2211 ierr = PetscLogObjectMemory((PetscObject)B,b->lda*b->Nmax*sizeof(PetscScalar));CHKERRQ(ierr); 2212 2213 b->user_alloc = PETSC_FALSE; 2214 } else { /* user-allocated storage */ 2215 if (!b->user_alloc) { ierr = PetscFree(b->v);CHKERRQ(ierr); } 2216 b->v = data; 2217 b->user_alloc = PETSC_TRUE; 2218 } 2219 B->assembled = PETSC_TRUE; 2220 PetscFunctionReturn(0); 2221 } 2222 2223 #undef __FUNCT__ 2224 #define __FUNCT__ "MatConvert_SeqDense_Elemental" 2225 PETSC_EXTERN PetscErrorCode MatConvert_SeqDense_Elemental(Mat A, MatType newtype,MatReuse reuse,Mat *newmat) 2226 { 2227 PetscFunctionBegin; 2228 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No support for requested operation"); 2229 PetscFunctionReturn(0); 2230 } 2231 2232 #undef __FUNCT__ 2233 #define __FUNCT__ "MatSeqDenseSetLDA" 2234 /*@C 2235 MatSeqDenseSetLDA - Declare the leading dimension of the user-provided array 2236 2237 Input parameter: 2238 + A - the matrix 2239 - lda - the leading dimension 2240 2241 Notes: 2242 This routine is to be used in conjunction with MatSeqDenseSetPreallocation(); 2243 it asserts that the preallocation has a leading dimension (the LDA parameter 2244 of Blas and Lapack fame) larger than M, the first dimension of the matrix. 2245 2246 Level: intermediate 2247 2248 .keywords: dense, matrix, LAPACK, BLAS 2249 2250 .seealso: MatCreate(), MatCreateSeqDense(), MatSeqDenseSetPreallocation(), MatSetMaximumSize() 2251 2252 @*/ 2253 PetscErrorCode MatSeqDenseSetLDA(Mat B,PetscInt lda) 2254 { 2255 Mat_SeqDense *b = (Mat_SeqDense*)B->data; 2256 2257 PetscFunctionBegin; 2258 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); 2259 b->lda = lda; 2260 b->changelda = PETSC_FALSE; 2261 b->Mmax = PetscMax(b->Mmax,lda); 2262 PetscFunctionReturn(0); 2263 } 2264 2265 /*MC 2266 MATSEQDENSE - MATSEQDENSE = "seqdense" - A matrix type to be used for sequential dense matrices. 2267 2268 Options Database Keys: 2269 . -mat_type seqdense - sets the matrix type to "seqdense" during a call to MatSetFromOptions() 2270 2271 Level: beginner 2272 2273 .seealso: MatCreateSeqDense() 2274 2275 M*/ 2276 2277 #undef __FUNCT__ 2278 #define __FUNCT__ "MatCreate_SeqDense" 2279 PETSC_EXTERN PetscErrorCode MatCreate_SeqDense(Mat B) 2280 { 2281 Mat_SeqDense *b; 2282 PetscErrorCode ierr; 2283 PetscMPIInt size; 2284 2285 PetscFunctionBegin; 2286 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)B),&size);CHKERRQ(ierr); 2287 if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Comm must be of size 1"); 2288 2289 ierr = PetscNewLog(B,&b);CHKERRQ(ierr); 2290 ierr = PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 2291 B->data = (void*)b; 2292 2293 b->pivots = 0; 2294 b->roworiented = PETSC_TRUE; 2295 b->v = 0; 2296 b->changelda = PETSC_FALSE; 2297 2298 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseGetArray_C",MatDenseGetArray_SeqDense);CHKERRQ(ierr); 2299 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreArray_C",MatDenseRestoreArray_SeqDense);CHKERRQ(ierr); 2300 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqdense_seqaij_C",MatConvert_SeqDense_SeqAIJ);CHKERRQ(ierr); 2301 #if defined(PETSC_HAVE_ELEMENTAL) 2302 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqdense_elemental_C",MatConvert_SeqDense_Elemental);CHKERRQ(ierr); 2303 #endif 2304 ierr = PetscObjectComposeFunction((PetscObject)B,"MatGetFactor_petsc_C",MatGetFactor_seqdense_petsc);CHKERRQ(ierr); 2305 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqDenseSetPreallocation_C",MatSeqDenseSetPreallocation_SeqDense);CHKERRQ(ierr); 2306 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMult_seqaij_seqdense_C",MatMatMult_SeqAIJ_SeqDense);CHKERRQ(ierr); 2307 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMultSymbolic_seqaij_seqdense_C",MatMatMultSymbolic_SeqAIJ_SeqDense);CHKERRQ(ierr); 2308 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMultNumeric_seqaij_seqdense_C",MatMatMultNumeric_SeqAIJ_SeqDense);CHKERRQ(ierr); 2309 ierr = PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMult_seqaij_seqdense_C",MatTransposeMatMult_SeqAIJ_SeqDense);CHKERRQ(ierr); 2310 ierr = PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultSymbolic_seqaij_seqdense_C",MatTransposeMatMultSymbolic_SeqAIJ_SeqDense);CHKERRQ(ierr); 2311 ierr = PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultNumeric_seqaij_seqdense_C",MatTransposeMatMultNumeric_SeqAIJ_SeqDense);CHKERRQ(ierr); 2312 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQDENSE);CHKERRQ(ierr); 2313 PetscFunctionReturn(0); 2314 } 2315