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