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