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