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