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