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