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