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