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