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