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