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 defined(PETSC_USE_COMPLEX) 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 #else 1068 if (v[j*m+i] > 0.) { 1069 color = PETSC_DRAW_RED; 1070 } else if (v[j*m+i] < 0.) { 1071 color = PETSC_DRAW_BLUE; 1072 } else { 1073 continue; 1074 } 1075 #endif 1076 ierr = PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);CHKERRQ(ierr); 1077 } 1078 } 1079 } else { 1080 /* use contour shading to indicate magnitude of values */ 1081 /* first determine max of all nonzero values */ 1082 for (i = 0; i < m*n; i++) { 1083 if (PetscAbsScalar(v[i]) > maxv) maxv = PetscAbsScalar(v[i]); 1084 } 1085 scale = (245.0 - PETSC_DRAW_BASIC_COLORS)/maxv; 1086 ierr = PetscDrawGetPopup(draw,&popup);CHKERRQ(ierr); 1087 if (popup) {ierr = PetscDrawScalePopup(popup,0.0,maxv);CHKERRQ(ierr);} 1088 for (j = 0; j < n; j++) { 1089 x_l = j; 1090 x_r = x_l + 1.0; 1091 for (i = 0; i < m; i++) { 1092 y_l = m - i - 1.0; 1093 y_r = y_l + 1.0; 1094 color = PETSC_DRAW_BASIC_COLORS + (int)(scale*PetscAbsScalar(v[j*m+i])); 1095 ierr = PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);CHKERRQ(ierr); 1096 } 1097 } 1098 } 1099 PetscFunctionReturn(0); 1100 } 1101 1102 #undef __FUNCT__ 1103 #define __FUNCT__ "MatView_SeqDense_Draw" 1104 PetscErrorCode MatView_SeqDense_Draw(Mat A,PetscViewer viewer) 1105 { 1106 PetscDraw draw; 1107 PetscBool isnull; 1108 PetscReal xr,yr,xl,yl,h,w; 1109 PetscErrorCode ierr; 1110 1111 PetscFunctionBegin; 1112 ierr = PetscViewerDrawGetDraw(viewer,0,&draw);CHKERRQ(ierr); 1113 ierr = PetscDrawIsNull(draw,&isnull);CHKERRQ(ierr); 1114 if (isnull) PetscFunctionReturn(0); 1115 1116 ierr = PetscObjectCompose((PetscObject)A,"Zoomviewer",(PetscObject)viewer);CHKERRQ(ierr); 1117 xr = A->cmap->n; yr = A->rmap->n; h = yr/10.0; w = xr/10.0; 1118 xr += w; yr += h; xl = -w; yl = -h; 1119 ierr = PetscDrawSetCoordinates(draw,xl,yl,xr,yr);CHKERRQ(ierr); 1120 ierr = PetscDrawZoom(draw,MatView_SeqDense_Draw_Zoom,A);CHKERRQ(ierr); 1121 ierr = PetscObjectCompose((PetscObject)A,"Zoomviewer",PETSC_NULL);CHKERRQ(ierr); 1122 PetscFunctionReturn(0); 1123 } 1124 1125 #undef __FUNCT__ 1126 #define __FUNCT__ "MatView_SeqDense" 1127 PetscErrorCode MatView_SeqDense(Mat A,PetscViewer viewer) 1128 { 1129 PetscErrorCode ierr; 1130 PetscBool iascii,isbinary,isdraw; 1131 1132 PetscFunctionBegin; 1133 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 1134 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);CHKERRQ(ierr); 1135 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);CHKERRQ(ierr); 1136 1137 if (iascii) { 1138 ierr = MatView_SeqDense_ASCII(A,viewer);CHKERRQ(ierr); 1139 } else if (isbinary) { 1140 ierr = MatView_SeqDense_Binary(A,viewer);CHKERRQ(ierr); 1141 } else if (isdraw) { 1142 ierr = MatView_SeqDense_Draw(A,viewer);CHKERRQ(ierr); 1143 } else { 1144 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Viewer type %s not supported by dense matrix",((PetscObject)viewer)->type_name); 1145 } 1146 PetscFunctionReturn(0); 1147 } 1148 1149 #undef __FUNCT__ 1150 #define __FUNCT__ "MatDestroy_SeqDense" 1151 PetscErrorCode MatDestroy_SeqDense(Mat mat) 1152 { 1153 Mat_SeqDense *l = (Mat_SeqDense*)mat->data; 1154 PetscErrorCode ierr; 1155 1156 PetscFunctionBegin; 1157 #if defined(PETSC_USE_LOG) 1158 PetscLogObjectState((PetscObject)mat,"Rows %D Cols %D",mat->rmap->n,mat->cmap->n); 1159 #endif 1160 ierr = PetscFree(l->pivots);CHKERRQ(ierr); 1161 if (!l->user_alloc) {ierr = PetscFree(l->v);CHKERRQ(ierr);} 1162 ierr = PetscFree(mat->data);CHKERRQ(ierr); 1163 1164 ierr = PetscObjectChangeTypeName((PetscObject)mat,0);CHKERRQ(ierr); 1165 ierr = PetscObjectComposeFunctionDynamic((PetscObject)mat,"MatDenseGetArray_C","",PETSC_NULL);CHKERRQ(ierr); 1166 ierr = PetscObjectComposeFunctionDynamic((PetscObject)mat,"MatDenseRestoreArray_C","",PETSC_NULL);CHKERRQ(ierr); 1167 ierr = PetscObjectComposeFunctionDynamic((PetscObject)mat,"MatSeqDenseSetPreallocation_C","",PETSC_NULL);CHKERRQ(ierr); 1168 ierr = PetscObjectComposeFunctionDynamic((PetscObject)mat,"MatMatMult_seqaij_seqdense_C","",PETSC_NULL);CHKERRQ(ierr); 1169 ierr = PetscObjectComposeFunctionDynamic((PetscObject)mat,"MatMatMultSymbolic_seqaij_seqdense_C","",PETSC_NULL);CHKERRQ(ierr); 1170 ierr = PetscObjectComposeFunctionDynamic((PetscObject)mat,"MatMatMultNumeric_seqaij_seqdense_C","",PETSC_NULL);CHKERRQ(ierr); 1171 PetscFunctionReturn(0); 1172 } 1173 1174 #undef __FUNCT__ 1175 #define __FUNCT__ "MatTranspose_SeqDense" 1176 PetscErrorCode MatTranspose_SeqDense(Mat A,MatReuse reuse,Mat *matout) 1177 { 1178 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1179 PetscErrorCode ierr; 1180 PetscInt k,j,m,n,M; 1181 PetscScalar *v,tmp; 1182 1183 PetscFunctionBegin; 1184 v = mat->v; m = A->rmap->n; M = mat->lda; n = A->cmap->n; 1185 if (reuse == MAT_REUSE_MATRIX && *matout == A) { /* in place transpose */ 1186 if (m != n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Can not transpose non-square matrix in place"); 1187 else { 1188 for (j=0; j<m; j++) { 1189 for (k=0; k<j; k++) { 1190 tmp = v[j + k*M]; 1191 v[j + k*M] = v[k + j*M]; 1192 v[k + j*M] = tmp; 1193 } 1194 } 1195 } 1196 } else { /* out-of-place transpose */ 1197 Mat tmat; 1198 Mat_SeqDense *tmatd; 1199 PetscScalar *v2; 1200 1201 if (reuse == MAT_INITIAL_MATRIX) { 1202 ierr = MatCreate(((PetscObject)A)->comm,&tmat);CHKERRQ(ierr); 1203 ierr = MatSetSizes(tmat,A->cmap->n,A->rmap->n,A->cmap->n,A->rmap->n);CHKERRQ(ierr); 1204 ierr = MatSetType(tmat,((PetscObject)A)->type_name);CHKERRQ(ierr); 1205 ierr = MatSeqDenseSetPreallocation(tmat,PETSC_NULL);CHKERRQ(ierr); 1206 } else { 1207 tmat = *matout; 1208 } 1209 tmatd = (Mat_SeqDense*)tmat->data; 1210 v = mat->v; v2 = tmatd->v; 1211 for (j=0; j<n; j++) { 1212 for (k=0; k<m; k++) v2[j + k*n] = v[k + j*M]; 1213 } 1214 ierr = MatAssemblyBegin(tmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1215 ierr = MatAssemblyEnd(tmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1216 *matout = tmat; 1217 } 1218 PetscFunctionReturn(0); 1219 } 1220 1221 #undef __FUNCT__ 1222 #define __FUNCT__ "MatEqual_SeqDense" 1223 PetscErrorCode MatEqual_SeqDense(Mat A1,Mat A2,PetscBool *flg) 1224 { 1225 Mat_SeqDense *mat1 = (Mat_SeqDense*)A1->data; 1226 Mat_SeqDense *mat2 = (Mat_SeqDense*)A2->data; 1227 PetscInt i,j; 1228 PetscScalar *v1,*v2; 1229 1230 PetscFunctionBegin; 1231 if (A1->rmap->n != A2->rmap->n) {*flg = PETSC_FALSE; PetscFunctionReturn(0);} 1232 if (A1->cmap->n != A2->cmap->n) {*flg = PETSC_FALSE; PetscFunctionReturn(0);} 1233 for (i=0; i<A1->rmap->n; i++) { 1234 v1 = mat1->v+i; v2 = mat2->v+i; 1235 for (j=0; j<A1->cmap->n; j++) { 1236 if (*v1 != *v2) {*flg = PETSC_FALSE; PetscFunctionReturn(0);} 1237 v1 += mat1->lda; v2 += mat2->lda; 1238 } 1239 } 1240 *flg = PETSC_TRUE; 1241 PetscFunctionReturn(0); 1242 } 1243 1244 #undef __FUNCT__ 1245 #define __FUNCT__ "MatGetDiagonal_SeqDense" 1246 PetscErrorCode MatGetDiagonal_SeqDense(Mat A,Vec v) 1247 { 1248 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1249 PetscErrorCode ierr; 1250 PetscInt i,n,len; 1251 PetscScalar *x,zero = 0.0; 1252 1253 PetscFunctionBegin; 1254 ierr = VecSet(v,zero);CHKERRQ(ierr); 1255 ierr = VecGetSize(v,&n);CHKERRQ(ierr); 1256 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 1257 len = PetscMin(A->rmap->n,A->cmap->n); 1258 if (n != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming mat and vec"); 1259 for (i=0; i<len; i++) { 1260 x[i] = mat->v[i*mat->lda + i]; 1261 } 1262 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 1263 PetscFunctionReturn(0); 1264 } 1265 1266 #undef __FUNCT__ 1267 #define __FUNCT__ "MatDiagonalScale_SeqDense" 1268 PetscErrorCode MatDiagonalScale_SeqDense(Mat A,Vec ll,Vec rr) 1269 { 1270 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1271 PetscScalar *l,*r,x,*v; 1272 PetscErrorCode ierr; 1273 PetscInt i,j,m = A->rmap->n,n = A->cmap->n; 1274 1275 PetscFunctionBegin; 1276 if (ll) { 1277 ierr = VecGetSize(ll,&m);CHKERRQ(ierr); 1278 ierr = VecGetArray(ll,&l);CHKERRQ(ierr); 1279 if (m != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Left scaling vec wrong size"); 1280 for (i=0; i<m; i++) { 1281 x = l[i]; 1282 v = mat->v + i; 1283 for (j=0; j<n; j++) { (*v) *= x; v+= m;} 1284 } 1285 ierr = VecRestoreArray(ll,&l);CHKERRQ(ierr); 1286 ierr = PetscLogFlops(n*m);CHKERRQ(ierr); 1287 } 1288 if (rr) { 1289 ierr = VecGetSize(rr,&n);CHKERRQ(ierr); 1290 ierr = VecGetArray(rr,&r);CHKERRQ(ierr); 1291 if (n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Right scaling vec wrong size"); 1292 for (i=0; i<n; i++) { 1293 x = r[i]; 1294 v = mat->v + i*m; 1295 for (j=0; j<m; j++) { (*v++) *= x;} 1296 } 1297 ierr = VecRestoreArray(rr,&r);CHKERRQ(ierr); 1298 ierr = PetscLogFlops(n*m);CHKERRQ(ierr); 1299 } 1300 PetscFunctionReturn(0); 1301 } 1302 1303 #undef __FUNCT__ 1304 #define __FUNCT__ "MatNorm_SeqDense" 1305 PetscErrorCode MatNorm_SeqDense(Mat A,NormType type,PetscReal *nrm) 1306 { 1307 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1308 PetscScalar *v = mat->v; 1309 PetscReal sum = 0.0; 1310 PetscInt lda=mat->lda,m=A->rmap->n,i,j; 1311 PetscErrorCode ierr; 1312 1313 PetscFunctionBegin; 1314 if (type == NORM_FROBENIUS) { 1315 if (lda>m) { 1316 for (j=0; j<A->cmap->n; j++) { 1317 v = mat->v+j*lda; 1318 for (i=0; i<m; i++) { 1319 sum += PetscRealPart(PetscConj(*v)*(*v)); v++; 1320 } 1321 } 1322 } else { 1323 for (i=0; i<A->cmap->n*A->rmap->n; i++) { 1324 sum += PetscRealPart(PetscConj(*v)*(*v)); v++; 1325 } 1326 } 1327 *nrm = PetscSqrtReal(sum); 1328 ierr = PetscLogFlops(2.0*A->cmap->n*A->rmap->n);CHKERRQ(ierr); 1329 } else if (type == NORM_1) { 1330 *nrm = 0.0; 1331 for (j=0; j<A->cmap->n; j++) { 1332 v = mat->v + j*mat->lda; 1333 sum = 0.0; 1334 for (i=0; i<A->rmap->n; i++) { 1335 sum += PetscAbsScalar(*v); v++; 1336 } 1337 if (sum > *nrm) *nrm = sum; 1338 } 1339 ierr = PetscLogFlops(A->cmap->n*A->rmap->n);CHKERRQ(ierr); 1340 } else if (type == NORM_INFINITY) { 1341 *nrm = 0.0; 1342 for (j=0; j<A->rmap->n; j++) { 1343 v = mat->v + j; 1344 sum = 0.0; 1345 for (i=0; i<A->cmap->n; i++) { 1346 sum += PetscAbsScalar(*v); v += mat->lda; 1347 } 1348 if (sum > *nrm) *nrm = sum; 1349 } 1350 ierr = PetscLogFlops(A->cmap->n*A->rmap->n);CHKERRQ(ierr); 1351 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No two norm"); 1352 PetscFunctionReturn(0); 1353 } 1354 1355 #undef __FUNCT__ 1356 #define __FUNCT__ "MatSetOption_SeqDense" 1357 PetscErrorCode MatSetOption_SeqDense(Mat A,MatOption op,PetscBool flg) 1358 { 1359 Mat_SeqDense *aij = (Mat_SeqDense*)A->data; 1360 PetscErrorCode ierr; 1361 1362 PetscFunctionBegin; 1363 switch (op) { 1364 case MAT_ROW_ORIENTED: 1365 aij->roworiented = flg; 1366 break; 1367 case MAT_NEW_NONZERO_LOCATIONS: 1368 case MAT_NEW_NONZERO_LOCATION_ERR: 1369 case MAT_NEW_NONZERO_ALLOCATION_ERR: 1370 case MAT_NEW_DIAGONALS: 1371 case MAT_KEEP_NONZERO_PATTERN: 1372 case MAT_IGNORE_OFF_PROC_ENTRIES: 1373 case MAT_USE_HASH_TABLE: 1374 case MAT_IGNORE_LOWER_TRIANGULAR: 1375 ierr = PetscInfo1(A,"Option %s ignored\n",MatOptions[op]);CHKERRQ(ierr); 1376 break; 1377 case MAT_SPD: 1378 case MAT_SYMMETRIC: 1379 case MAT_STRUCTURALLY_SYMMETRIC: 1380 case MAT_HERMITIAN: 1381 case MAT_SYMMETRY_ETERNAL: 1382 /* These options are handled directly by MatSetOption() */ 1383 break; 1384 default: 1385 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"unknown option %s",MatOptions[op]); 1386 } 1387 PetscFunctionReturn(0); 1388 } 1389 1390 #undef __FUNCT__ 1391 #define __FUNCT__ "MatZeroEntries_SeqDense" 1392 PetscErrorCode MatZeroEntries_SeqDense(Mat A) 1393 { 1394 Mat_SeqDense *l = (Mat_SeqDense*)A->data; 1395 PetscErrorCode ierr; 1396 PetscInt lda=l->lda,m=A->rmap->n,j; 1397 1398 PetscFunctionBegin; 1399 if (lda>m) { 1400 for (j=0; j<A->cmap->n; j++) { 1401 ierr = PetscMemzero(l->v+j*lda,m*sizeof(PetscScalar));CHKERRQ(ierr); 1402 } 1403 } else { 1404 ierr = PetscMemzero(l->v,A->rmap->n*A->cmap->n*sizeof(PetscScalar));CHKERRQ(ierr); 1405 } 1406 PetscFunctionReturn(0); 1407 } 1408 1409 #undef __FUNCT__ 1410 #define __FUNCT__ "MatZeroRows_SeqDense" 1411 PetscErrorCode MatZeroRows_SeqDense(Mat A,PetscInt N,const PetscInt rows[],PetscScalar diag,Vec x,Vec b) 1412 { 1413 PetscErrorCode ierr; 1414 Mat_SeqDense *l = (Mat_SeqDense*)A->data; 1415 PetscInt m = l->lda, n = A->cmap->n, i,j; 1416 PetscScalar *slot,*bb; 1417 const PetscScalar *xx; 1418 1419 PetscFunctionBegin; 1420 #if defined(PETSC_USE_DEBUG) 1421 for (i=0; i<N; i++) { 1422 if (rows[i] < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative row requested to be zeroed"); 1423 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); 1424 } 1425 #endif 1426 1427 /* fix right hand side if needed */ 1428 if (x && b) { 1429 ierr = VecGetArrayRead(x,&xx);CHKERRQ(ierr); 1430 ierr = VecGetArray(b,&bb);CHKERRQ(ierr); 1431 for (i=0; i<N; i++) { 1432 bb[rows[i]] = diag*xx[rows[i]]; 1433 } 1434 ierr = VecRestoreArrayRead(x,&xx);CHKERRQ(ierr); 1435 ierr = VecRestoreArray(b,&bb);CHKERRQ(ierr); 1436 } 1437 1438 for (i=0; i<N; i++) { 1439 slot = l->v + rows[i]; 1440 for (j=0; j<n; j++) { *slot = 0.0; slot += m;} 1441 } 1442 if (diag != 0.0) { 1443 if (A->rmap->n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only coded for square matrices"); 1444 for (i=0; i<N; i++) { 1445 slot = l->v + (m+1)*rows[i]; 1446 *slot = diag; 1447 } 1448 } 1449 PetscFunctionReturn(0); 1450 } 1451 1452 #undef __FUNCT__ 1453 #define __FUNCT__ "MatDenseGetArray_SeqDense" 1454 PetscErrorCode MatDenseGetArray_SeqDense(Mat A,PetscScalar *array[]) 1455 { 1456 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1457 1458 PetscFunctionBegin; 1459 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"); 1460 *array = mat->v; 1461 PetscFunctionReturn(0); 1462 } 1463 1464 #undef __FUNCT__ 1465 #define __FUNCT__ "MatDenseRestoreArray_SeqDense" 1466 PetscErrorCode MatDenseRestoreArray_SeqDense(Mat A,PetscScalar *array[]) 1467 { 1468 PetscFunctionBegin; 1469 *array = 0; /* user cannot accidently use the array later */ 1470 PetscFunctionReturn(0); 1471 } 1472 1473 #undef __FUNCT__ 1474 #define __FUNCT__ "MatDenseGetArray" 1475 /*@C 1476 MatDenseGetArray - gives access to the array where the data for a SeqDense matrix is stored 1477 1478 Not Collective 1479 1480 Input Parameter: 1481 . mat - a MATSEQDENSE matrix 1482 1483 Output Parameter: 1484 . array - pointer to the data 1485 1486 Level: intermediate 1487 1488 .seealso: MatDenseRestoreArray() 1489 @*/ 1490 PetscErrorCode MatDenseGetArray(Mat A,PetscScalar **array) 1491 { 1492 PetscErrorCode ierr; 1493 1494 PetscFunctionBegin; 1495 ierr = PetscUseMethod(A,"MatDenseGetArray_C",(Mat,PetscScalar**),(A,array));CHKERRQ(ierr); 1496 PetscFunctionReturn(0); 1497 } 1498 1499 #undef __FUNCT__ 1500 #define __FUNCT__ "MatDenseRestoreArray" 1501 /*@C 1502 MatDenseRestoreArray - returns access to the array where the data for a SeqDense matrix is stored obtained by MatDenseGetArray() 1503 1504 Not Collective 1505 1506 Input Parameters: 1507 . mat - a MATSEQDENSE matrix 1508 . array - pointer to the data 1509 1510 Level: intermediate 1511 1512 .seealso: MatDenseGetArray() 1513 @*/ 1514 PetscErrorCode MatDenseRestoreArray(Mat A,PetscScalar **array) 1515 { 1516 PetscErrorCode ierr; 1517 1518 PetscFunctionBegin; 1519 ierr = PetscUseMethod(A,"MatDenseRestoreArray_C",(Mat,PetscScalar**),(A,array));CHKERRQ(ierr); 1520 PetscFunctionReturn(0); 1521 } 1522 1523 #undef __FUNCT__ 1524 #define __FUNCT__ "MatGetSubMatrix_SeqDense" 1525 static PetscErrorCode MatGetSubMatrix_SeqDense(Mat A,IS isrow,IS iscol,PetscInt cs,MatReuse scall,Mat *B) 1526 { 1527 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1528 PetscErrorCode ierr; 1529 PetscInt i,j,nrows,ncols; 1530 const PetscInt *irow,*icol; 1531 PetscScalar *av,*bv,*v = mat->v; 1532 Mat newmat; 1533 1534 PetscFunctionBegin; 1535 ierr = ISGetIndices(isrow,&irow);CHKERRQ(ierr); 1536 ierr = ISGetIndices(iscol,&icol);CHKERRQ(ierr); 1537 ierr = ISGetLocalSize(isrow,&nrows);CHKERRQ(ierr); 1538 ierr = ISGetLocalSize(iscol,&ncols);CHKERRQ(ierr); 1539 1540 /* Check submatrixcall */ 1541 if (scall == MAT_REUSE_MATRIX) { 1542 PetscInt n_cols,n_rows; 1543 ierr = MatGetSize(*B,&n_rows,&n_cols);CHKERRQ(ierr); 1544 if (n_rows != nrows || n_cols != ncols) { 1545 /* resize the result matrix to match number of requested rows/columns */ 1546 ierr = MatSetSizes(*B,nrows,ncols,nrows,ncols);CHKERRQ(ierr); 1547 } 1548 newmat = *B; 1549 } else { 1550 /* Create and fill new matrix */ 1551 ierr = MatCreate(((PetscObject)A)->comm,&newmat);CHKERRQ(ierr); 1552 ierr = MatSetSizes(newmat,nrows,ncols,nrows,ncols);CHKERRQ(ierr); 1553 ierr = MatSetType(newmat,((PetscObject)A)->type_name);CHKERRQ(ierr); 1554 ierr = MatSeqDenseSetPreallocation(newmat,PETSC_NULL);CHKERRQ(ierr); 1555 } 1556 1557 /* Now extract the data pointers and do the copy,column at a time */ 1558 bv = ((Mat_SeqDense*)newmat->data)->v; 1559 1560 for (i=0; i<ncols; i++) { 1561 av = v + mat->lda*icol[i]; 1562 for (j=0; j<nrows; j++) { 1563 *bv++ = av[irow[j]]; 1564 } 1565 } 1566 1567 /* Assemble the matrices so that the correct flags are set */ 1568 ierr = MatAssemblyBegin(newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1569 ierr = MatAssemblyEnd(newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1570 1571 /* Free work space */ 1572 ierr = ISRestoreIndices(isrow,&irow);CHKERRQ(ierr); 1573 ierr = ISRestoreIndices(iscol,&icol);CHKERRQ(ierr); 1574 *B = newmat; 1575 PetscFunctionReturn(0); 1576 } 1577 1578 #undef __FUNCT__ 1579 #define __FUNCT__ "MatGetSubMatrices_SeqDense" 1580 PetscErrorCode MatGetSubMatrices_SeqDense(Mat A,PetscInt n,const IS irow[],const IS icol[],MatReuse scall,Mat *B[]) 1581 { 1582 PetscErrorCode ierr; 1583 PetscInt i; 1584 1585 PetscFunctionBegin; 1586 if (scall == MAT_INITIAL_MATRIX) { 1587 ierr = PetscMalloc((n+1)*sizeof(Mat),B);CHKERRQ(ierr); 1588 } 1589 1590 for (i=0; i<n; i++) { 1591 ierr = MatGetSubMatrix_SeqDense(A,irow[i],icol[i],PETSC_DECIDE,scall,&(*B)[i]);CHKERRQ(ierr); 1592 } 1593 PetscFunctionReturn(0); 1594 } 1595 1596 #undef __FUNCT__ 1597 #define __FUNCT__ "MatAssemblyBegin_SeqDense" 1598 PetscErrorCode MatAssemblyBegin_SeqDense(Mat mat,MatAssemblyType mode) 1599 { 1600 PetscFunctionBegin; 1601 PetscFunctionReturn(0); 1602 } 1603 1604 #undef __FUNCT__ 1605 #define __FUNCT__ "MatAssemblyEnd_SeqDense" 1606 PetscErrorCode MatAssemblyEnd_SeqDense(Mat mat,MatAssemblyType mode) 1607 { 1608 PetscFunctionBegin; 1609 PetscFunctionReturn(0); 1610 } 1611 1612 #undef __FUNCT__ 1613 #define __FUNCT__ "MatCopy_SeqDense" 1614 PetscErrorCode MatCopy_SeqDense(Mat A,Mat B,MatStructure str) 1615 { 1616 Mat_SeqDense *a = (Mat_SeqDense*)A->data,*b = (Mat_SeqDense *)B->data; 1617 PetscErrorCode ierr; 1618 PetscInt lda1=a->lda,lda2=b->lda, m=A->rmap->n,n=A->cmap->n, j; 1619 1620 PetscFunctionBegin; 1621 /* If the two matrices don't have the same copy implementation, they aren't compatible for fast copy. */ 1622 if (A->ops->copy != B->ops->copy) { 1623 ierr = MatCopy_Basic(A,B,str);CHKERRQ(ierr); 1624 PetscFunctionReturn(0); 1625 } 1626 if (m != B->rmap->n || n != B->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"size(B) != size(A)"); 1627 if (lda1>m || lda2>m) { 1628 for (j=0; j<n; j++) { 1629 ierr = PetscMemcpy(b->v+j*lda2,a->v+j*lda1,m*sizeof(PetscScalar));CHKERRQ(ierr); 1630 } 1631 } else { 1632 ierr = PetscMemcpy(b->v,a->v,A->rmap->n*A->cmap->n*sizeof(PetscScalar));CHKERRQ(ierr); 1633 } 1634 PetscFunctionReturn(0); 1635 } 1636 1637 #undef __FUNCT__ 1638 #define __FUNCT__ "MatSetUp_SeqDense" 1639 PetscErrorCode MatSetUp_SeqDense(Mat A) 1640 { 1641 PetscErrorCode ierr; 1642 1643 PetscFunctionBegin; 1644 ierr = MatSeqDenseSetPreallocation(A,0);CHKERRQ(ierr); 1645 PetscFunctionReturn(0); 1646 } 1647 1648 #undef __FUNCT__ 1649 #define __FUNCT__ "MatConjugate_SeqDense" 1650 static PetscErrorCode MatConjugate_SeqDense(Mat A) 1651 { 1652 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1653 PetscInt i,nz = A->rmap->n*A->cmap->n; 1654 PetscScalar *aa = a->v; 1655 1656 PetscFunctionBegin; 1657 for (i=0; i<nz; i++) aa[i] = PetscConj(aa[i]); 1658 PetscFunctionReturn(0); 1659 } 1660 1661 #undef __FUNCT__ 1662 #define __FUNCT__ "MatRealPart_SeqDense" 1663 static PetscErrorCode MatRealPart_SeqDense(Mat A) 1664 { 1665 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1666 PetscInt i,nz = A->rmap->n*A->cmap->n; 1667 PetscScalar *aa = a->v; 1668 1669 PetscFunctionBegin; 1670 for (i=0; i<nz; i++) aa[i] = PetscRealPart(aa[i]); 1671 PetscFunctionReturn(0); 1672 } 1673 1674 #undef __FUNCT__ 1675 #define __FUNCT__ "MatImaginaryPart_SeqDense" 1676 static PetscErrorCode MatImaginaryPart_SeqDense(Mat A) 1677 { 1678 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1679 PetscInt i,nz = A->rmap->n*A->cmap->n; 1680 PetscScalar *aa = a->v; 1681 1682 PetscFunctionBegin; 1683 for (i=0; i<nz; i++) aa[i] = PetscImaginaryPart(aa[i]); 1684 PetscFunctionReturn(0); 1685 } 1686 1687 /* ----------------------------------------------------------------*/ 1688 #undef __FUNCT__ 1689 #define __FUNCT__ "MatMatMult_SeqDense_SeqDense" 1690 PetscErrorCode MatMatMult_SeqDense_SeqDense(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C) 1691 { 1692 PetscErrorCode ierr; 1693 1694 PetscFunctionBegin; 1695 if (scall == MAT_INITIAL_MATRIX){ 1696 ierr = MatMatMultSymbolic_SeqDense_SeqDense(A,B,fill,C);CHKERRQ(ierr); 1697 } 1698 ierr = MatMatMultNumeric_SeqDense_SeqDense(A,B,*C);CHKERRQ(ierr); 1699 PetscFunctionReturn(0); 1700 } 1701 1702 #undef __FUNCT__ 1703 #define __FUNCT__ "MatMatMultSymbolic_SeqDense_SeqDense" 1704 PetscErrorCode MatMatMultSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat *C) 1705 { 1706 PetscErrorCode ierr; 1707 PetscInt m=A->rmap->n,n=B->cmap->n; 1708 Mat Cmat; 1709 1710 PetscFunctionBegin; 1711 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); 1712 ierr = MatCreate(PETSC_COMM_SELF,&Cmat);CHKERRQ(ierr); 1713 ierr = MatSetSizes(Cmat,m,n,m,n);CHKERRQ(ierr); 1714 ierr = MatSetType(Cmat,MATSEQDENSE);CHKERRQ(ierr); 1715 ierr = MatSeqDenseSetPreallocation(Cmat,PETSC_NULL);CHKERRQ(ierr); 1716 1717 *C = Cmat; 1718 PetscFunctionReturn(0); 1719 } 1720 1721 #undef __FUNCT__ 1722 #define __FUNCT__ "MatMatMultNumeric_SeqDense_SeqDense" 1723 PetscErrorCode MatMatMultNumeric_SeqDense_SeqDense(Mat A,Mat B,Mat C) 1724 { 1725 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1726 Mat_SeqDense *b = (Mat_SeqDense*)B->data; 1727 Mat_SeqDense *c = (Mat_SeqDense*)C->data; 1728 PetscBLASInt m,n,k; 1729 PetscScalar _DOne=1.0,_DZero=0.0; 1730 1731 PetscFunctionBegin; 1732 m = PetscBLASIntCast(A->rmap->n); 1733 n = PetscBLASIntCast(B->cmap->n); 1734 k = PetscBLASIntCast(A->cmap->n); 1735 BLASgemm_("N","N",&m,&n,&k,&_DOne,a->v,&a->lda,b->v,&b->lda,&_DZero,c->v,&c->lda); 1736 PetscFunctionReturn(0); 1737 } 1738 1739 #undef __FUNCT__ 1740 #define __FUNCT__ "MatTransposeMatMult_SeqDense_SeqDense" 1741 PetscErrorCode MatTransposeMatMult_SeqDense_SeqDense(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C) 1742 { 1743 PetscErrorCode ierr; 1744 1745 PetscFunctionBegin; 1746 if (scall == MAT_INITIAL_MATRIX){ 1747 ierr = MatTransposeMatMultSymbolic_SeqDense_SeqDense(A,B,fill,C);CHKERRQ(ierr); 1748 } 1749 ierr = MatTransposeMatMultNumeric_SeqDense_SeqDense(A,B,*C);CHKERRQ(ierr); 1750 PetscFunctionReturn(0); 1751 } 1752 1753 #undef __FUNCT__ 1754 #define __FUNCT__ "MatTransposeMatMultSymbolic_SeqDense_SeqDense" 1755 PetscErrorCode MatTransposeMatMultSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat *C) 1756 { 1757 PetscErrorCode ierr; 1758 PetscInt m=A->cmap->n,n=B->cmap->n; 1759 Mat Cmat; 1760 1761 PetscFunctionBegin; 1762 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); 1763 ierr = MatCreate(PETSC_COMM_SELF,&Cmat);CHKERRQ(ierr); 1764 ierr = MatSetSizes(Cmat,m,n,m,n);CHKERRQ(ierr); 1765 ierr = MatSetType(Cmat,MATSEQDENSE);CHKERRQ(ierr); 1766 ierr = MatSeqDenseSetPreallocation(Cmat,PETSC_NULL);CHKERRQ(ierr); 1767 Cmat->assembled = PETSC_TRUE; 1768 *C = Cmat; 1769 PetscFunctionReturn(0); 1770 } 1771 1772 #undef __FUNCT__ 1773 #define __FUNCT__ "MatTransposeMatMultNumeric_SeqDense_SeqDense" 1774 PetscErrorCode MatTransposeMatMultNumeric_SeqDense_SeqDense(Mat A,Mat B,Mat C) 1775 { 1776 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1777 Mat_SeqDense *b = (Mat_SeqDense*)B->data; 1778 Mat_SeqDense *c = (Mat_SeqDense*)C->data; 1779 PetscBLASInt m,n,k; 1780 PetscScalar _DOne=1.0,_DZero=0.0; 1781 1782 PetscFunctionBegin; 1783 m = PetscBLASIntCast(A->cmap->n); 1784 n = PetscBLASIntCast(B->cmap->n); 1785 k = PetscBLASIntCast(A->rmap->n); 1786 /* 1787 Note the m and n arguments below are the number rows and columns of A', not A! 1788 */ 1789 BLASgemm_("T","N",&m,&n,&k,&_DOne,a->v,&a->lda,b->v,&b->lda,&_DZero,c->v,&c->lda); 1790 PetscFunctionReturn(0); 1791 } 1792 1793 #undef __FUNCT__ 1794 #define __FUNCT__ "MatGetRowMax_SeqDense" 1795 PetscErrorCode MatGetRowMax_SeqDense(Mat A,Vec v,PetscInt idx[]) 1796 { 1797 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1798 PetscErrorCode ierr; 1799 PetscInt i,j,m = A->rmap->n,n = A->cmap->n,p; 1800 PetscScalar *x; 1801 MatScalar *aa = a->v; 1802 1803 PetscFunctionBegin; 1804 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 1805 1806 ierr = VecSet(v,0.0);CHKERRQ(ierr); 1807 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 1808 ierr = VecGetLocalSize(v,&p);CHKERRQ(ierr); 1809 if (p != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 1810 for (i=0; i<m; i++) { 1811 x[i] = aa[i]; if (idx) idx[i] = 0; 1812 for (j=1; j<n; j++){ 1813 if (PetscRealPart(x[i]) < PetscRealPart(aa[i+m*j])) {x[i] = aa[i + m*j]; if (idx) idx[i] = j;} 1814 } 1815 } 1816 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 1817 PetscFunctionReturn(0); 1818 } 1819 1820 #undef __FUNCT__ 1821 #define __FUNCT__ "MatGetRowMaxAbs_SeqDense" 1822 PetscErrorCode MatGetRowMaxAbs_SeqDense(Mat A,Vec v,PetscInt idx[]) 1823 { 1824 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1825 PetscErrorCode ierr; 1826 PetscInt i,j,m = A->rmap->n,n = A->cmap->n,p; 1827 PetscScalar *x; 1828 PetscReal atmp; 1829 MatScalar *aa = a->v; 1830 1831 PetscFunctionBegin; 1832 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 1833 1834 ierr = VecSet(v,0.0);CHKERRQ(ierr); 1835 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 1836 ierr = VecGetLocalSize(v,&p);CHKERRQ(ierr); 1837 if (p != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 1838 for (i=0; i<m; i++) { 1839 x[i] = PetscAbsScalar(aa[i]); 1840 for (j=1; j<n; j++){ 1841 atmp = PetscAbsScalar(aa[i+m*j]); 1842 if (PetscAbsScalar(x[i]) < atmp) {x[i] = atmp; if (idx) idx[i] = j;} 1843 } 1844 } 1845 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 1846 PetscFunctionReturn(0); 1847 } 1848 1849 #undef __FUNCT__ 1850 #define __FUNCT__ "MatGetRowMin_SeqDense" 1851 PetscErrorCode MatGetRowMin_SeqDense(Mat A,Vec v,PetscInt idx[]) 1852 { 1853 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1854 PetscErrorCode ierr; 1855 PetscInt i,j,m = A->rmap->n,n = A->cmap->n,p; 1856 PetscScalar *x; 1857 MatScalar *aa = a->v; 1858 1859 PetscFunctionBegin; 1860 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 1861 1862 ierr = VecSet(v,0.0);CHKERRQ(ierr); 1863 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 1864 ierr = VecGetLocalSize(v,&p);CHKERRQ(ierr); 1865 if (p != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 1866 for (i=0; i<m; i++) { 1867 x[i] = aa[i]; if (idx) idx[i] = 0; 1868 for (j=1; j<n; j++){ 1869 if (PetscRealPart(x[i]) > PetscRealPart(aa[i+m*j])) {x[i] = aa[i + m*j]; if (idx) idx[i] = j;} 1870 } 1871 } 1872 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 1873 PetscFunctionReturn(0); 1874 } 1875 1876 #undef __FUNCT__ 1877 #define __FUNCT__ "MatGetColumnVector_SeqDense" 1878 PetscErrorCode MatGetColumnVector_SeqDense(Mat A,Vec v,PetscInt col) 1879 { 1880 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1881 PetscErrorCode ierr; 1882 PetscScalar *x; 1883 1884 PetscFunctionBegin; 1885 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 1886 1887 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 1888 ierr = PetscMemcpy(x,a->v+col*a->lda,A->rmap->n*sizeof(PetscScalar));CHKERRQ(ierr); 1889 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 1890 PetscFunctionReturn(0); 1891 } 1892 1893 1894 #undef __FUNCT__ 1895 #define __FUNCT__ "MatGetColumnNorms_SeqDense" 1896 PetscErrorCode MatGetColumnNorms_SeqDense(Mat A,NormType type,PetscReal *norms) 1897 { 1898 PetscErrorCode ierr; 1899 PetscInt i,j,m,n; 1900 PetscScalar *a; 1901 1902 PetscFunctionBegin; 1903 ierr = MatGetSize(A,&m,&n);CHKERRQ(ierr); 1904 ierr = PetscMemzero(norms,n*sizeof(PetscReal));CHKERRQ(ierr); 1905 ierr = MatDenseGetArray(A,&a);CHKERRQ(ierr); 1906 if (type == NORM_2) { 1907 for (i=0; i<n; i++ ){ 1908 for (j=0; j<m; j++) { 1909 norms[i] += PetscAbsScalar(a[j]*a[j]); 1910 } 1911 a += m; 1912 } 1913 } else if (type == NORM_1) { 1914 for (i=0; i<n; i++ ){ 1915 for (j=0; j<m; j++) { 1916 norms[i] += PetscAbsScalar(a[j]); 1917 } 1918 a += m; 1919 } 1920 } else if (type == NORM_INFINITY) { 1921 for (i=0; i<n; i++ ){ 1922 for (j=0; j<m; j++) { 1923 norms[i] = PetscMax(PetscAbsScalar(a[j]),norms[i]); 1924 } 1925 a += m; 1926 } 1927 } else SETERRQ(((PetscObject)A)->comm,PETSC_ERR_ARG_WRONG,"Unknown NormType"); 1928 ierr = MatDenseRestoreArray(A,&a);CHKERRQ(ierr); 1929 if (type == NORM_2) { 1930 for (i=0; i<n; i++) norms[i] = PetscSqrtReal(norms[i]); 1931 } 1932 PetscFunctionReturn(0); 1933 } 1934 1935 #undef __FUNCT__ 1936 #define __FUNCT__ "MatSetRandom_SeqDense" 1937 static PetscErrorCode MatSetRandom_SeqDense(Mat x,PetscRandom rctx) 1938 { 1939 PetscErrorCode ierr; 1940 PetscScalar *a; 1941 PetscInt m,n,i; 1942 1943 PetscFunctionBegin; 1944 ierr = MatGetSize(x,&m,&n);CHKERRQ(ierr); 1945 ierr = MatDenseGetArray(x,&a);CHKERRQ(ierr); 1946 for (i=0; i<m*n; i++) { 1947 ierr = PetscRandomGetValue(rctx,a+i);CHKERRQ(ierr); 1948 } 1949 ierr = MatDenseRestoreArray(x,&a);CHKERRQ(ierr); 1950 PetscFunctionReturn(0); 1951 } 1952 1953 1954 /* -------------------------------------------------------------------*/ 1955 static struct _MatOps MatOps_Values = {MatSetValues_SeqDense, 1956 MatGetRow_SeqDense, 1957 MatRestoreRow_SeqDense, 1958 MatMult_SeqDense, 1959 /* 4*/ MatMultAdd_SeqDense, 1960 MatMultTranspose_SeqDense, 1961 MatMultTransposeAdd_SeqDense, 1962 0, 1963 0, 1964 0, 1965 /*10*/ 0, 1966 MatLUFactor_SeqDense, 1967 MatCholeskyFactor_SeqDense, 1968 MatSOR_SeqDense, 1969 MatTranspose_SeqDense, 1970 /*15*/ MatGetInfo_SeqDense, 1971 MatEqual_SeqDense, 1972 MatGetDiagonal_SeqDense, 1973 MatDiagonalScale_SeqDense, 1974 MatNorm_SeqDense, 1975 /*20*/ MatAssemblyBegin_SeqDense, 1976 MatAssemblyEnd_SeqDense, 1977 MatSetOption_SeqDense, 1978 MatZeroEntries_SeqDense, 1979 /*24*/ MatZeroRows_SeqDense, 1980 0, 1981 0, 1982 0, 1983 0, 1984 /*29*/ MatSetUp_SeqDense, 1985 0, 1986 0, 1987 0, 1988 0, 1989 /*34*/ MatDuplicate_SeqDense, 1990 0, 1991 0, 1992 0, 1993 0, 1994 /*39*/ MatAXPY_SeqDense, 1995 MatGetSubMatrices_SeqDense, 1996 0, 1997 MatGetValues_SeqDense, 1998 MatCopy_SeqDense, 1999 /*44*/ MatGetRowMax_SeqDense, 2000 MatScale_SeqDense, 2001 0, 2002 0, 2003 0, 2004 /*49*/ MatSetRandom_SeqDense, 2005 0, 2006 0, 2007 0, 2008 0, 2009 /*54*/ 0, 2010 0, 2011 0, 2012 0, 2013 0, 2014 /*59*/ 0, 2015 MatDestroy_SeqDense, 2016 MatView_SeqDense, 2017 0, 2018 0, 2019 /*64*/ 0, 2020 0, 2021 0, 2022 0, 2023 0, 2024 /*69*/ MatGetRowMaxAbs_SeqDense, 2025 0, 2026 0, 2027 0, 2028 0, 2029 /*74*/ 0, 2030 0, 2031 0, 2032 0, 2033 0, 2034 /*79*/ 0, 2035 0, 2036 0, 2037 0, 2038 /*83*/ MatLoad_SeqDense, 2039 0, 2040 MatIsHermitian_SeqDense, 2041 0, 2042 0, 2043 0, 2044 /*89*/ MatMatMult_SeqDense_SeqDense, 2045 MatMatMultSymbolic_SeqDense_SeqDense, 2046 MatMatMultNumeric_SeqDense_SeqDense, 2047 0, 2048 0, 2049 /*94*/ 0, 2050 0, 2051 0, 2052 0, 2053 0, 2054 /*99*/ 0, 2055 0, 2056 0, 2057 MatConjugate_SeqDense, 2058 0, 2059 /*104*/0, 2060 MatRealPart_SeqDense, 2061 MatImaginaryPart_SeqDense, 2062 0, 2063 0, 2064 /*109*/MatMatSolve_SeqDense, 2065 0, 2066 MatGetRowMin_SeqDense, 2067 MatGetColumnVector_SeqDense, 2068 0, 2069 /*114*/0, 2070 0, 2071 0, 2072 0, 2073 0, 2074 /*119*/0, 2075 0, 2076 0, 2077 0, 2078 0, 2079 /*124*/0, 2080 MatGetColumnNorms_SeqDense, 2081 0, 2082 0, 2083 0, 2084 /*129*/0, 2085 MatTransposeMatMult_SeqDense_SeqDense, 2086 MatTransposeMatMultSymbolic_SeqDense_SeqDense, 2087 MatTransposeMatMultNumeric_SeqDense_SeqDense, 2088 }; 2089 2090 #undef __FUNCT__ 2091 #define __FUNCT__ "MatCreateSeqDense" 2092 /*@C 2093 MatCreateSeqDense - Creates a sequential dense matrix that 2094 is stored in column major order (the usual Fortran 77 manner). Many 2095 of the matrix operations use the BLAS and LAPACK routines. 2096 2097 Collective on MPI_Comm 2098 2099 Input Parameters: 2100 + comm - MPI communicator, set to PETSC_COMM_SELF 2101 . m - number of rows 2102 . n - number of columns 2103 - data - optional location of matrix data in column major order. Set data=PETSC_NULL for PETSc 2104 to control all matrix memory allocation. 2105 2106 Output Parameter: 2107 . A - the matrix 2108 2109 Notes: 2110 The data input variable is intended primarily for Fortran programmers 2111 who wish to allocate their own matrix memory space. Most users should 2112 set data=PETSC_NULL. 2113 2114 Level: intermediate 2115 2116 .keywords: dense, matrix, LAPACK, BLAS 2117 2118 .seealso: MatCreate(), MatCreateDense(), MatSetValues() 2119 @*/ 2120 PetscErrorCode MatCreateSeqDense(MPI_Comm comm,PetscInt m,PetscInt n,PetscScalar *data,Mat *A) 2121 { 2122 PetscErrorCode ierr; 2123 2124 PetscFunctionBegin; 2125 ierr = MatCreate(comm,A);CHKERRQ(ierr); 2126 ierr = MatSetSizes(*A,m,n,m,n);CHKERRQ(ierr); 2127 ierr = MatSetType(*A,MATSEQDENSE);CHKERRQ(ierr); 2128 ierr = MatSeqDenseSetPreallocation(*A,data);CHKERRQ(ierr); 2129 PetscFunctionReturn(0); 2130 } 2131 2132 #undef __FUNCT__ 2133 #define __FUNCT__ "MatSeqDenseSetPreallocation" 2134 /*@C 2135 MatSeqDenseSetPreallocation - Sets the array used for storing the matrix elements 2136 2137 Collective on MPI_Comm 2138 2139 Input Parameters: 2140 + A - the matrix 2141 - data - the array (or PETSC_NULL) 2142 2143 Notes: 2144 The data input variable is intended primarily for Fortran programmers 2145 who wish to allocate their own matrix memory space. Most users should 2146 need not call this routine. 2147 2148 Level: intermediate 2149 2150 .keywords: dense, matrix, LAPACK, BLAS 2151 2152 .seealso: MatCreate(), MatCreateDense(), MatSetValues(), MatSeqDenseSetLDA() 2153 2154 @*/ 2155 PetscErrorCode MatSeqDenseSetPreallocation(Mat B,PetscScalar data[]) 2156 { 2157 PetscErrorCode ierr; 2158 2159 PetscFunctionBegin; 2160 ierr = PetscTryMethod(B,"MatSeqDenseSetPreallocation_C",(Mat,PetscScalar[]),(B,data));CHKERRQ(ierr); 2161 PetscFunctionReturn(0); 2162 } 2163 2164 EXTERN_C_BEGIN 2165 #undef __FUNCT__ 2166 #define __FUNCT__ "MatSeqDenseSetPreallocation_SeqDense" 2167 PetscErrorCode MatSeqDenseSetPreallocation_SeqDense(Mat B,PetscScalar *data) 2168 { 2169 Mat_SeqDense *b; 2170 PetscErrorCode ierr; 2171 2172 PetscFunctionBegin; 2173 B->preallocated = PETSC_TRUE; 2174 2175 ierr = PetscLayoutSetUp(B->rmap);CHKERRQ(ierr); 2176 ierr = PetscLayoutSetUp(B->cmap);CHKERRQ(ierr); 2177 2178 b = (Mat_SeqDense*)B->data; 2179 b->Mmax = B->rmap->n; 2180 b->Nmax = B->cmap->n; 2181 if (b->lda <= 0 || b->changelda) b->lda = B->rmap->n; 2182 2183 if (!data) { /* petsc-allocated storage */ 2184 if (!b->user_alloc) { ierr = PetscFree(b->v);CHKERRQ(ierr); } 2185 ierr = PetscMalloc(b->lda*b->Nmax*sizeof(PetscScalar),&b->v);CHKERRQ(ierr); 2186 ierr = PetscMemzero(b->v,b->lda*b->Nmax*sizeof(PetscScalar));CHKERRQ(ierr); 2187 ierr = PetscLogObjectMemory(B,b->lda*b->Nmax*sizeof(PetscScalar));CHKERRQ(ierr); 2188 b->user_alloc = PETSC_FALSE; 2189 } else { /* user-allocated storage */ 2190 if (!b->user_alloc) { ierr = PetscFree(b->v);CHKERRQ(ierr); } 2191 b->v = data; 2192 b->user_alloc = PETSC_TRUE; 2193 } 2194 B->assembled = PETSC_TRUE; 2195 PetscFunctionReturn(0); 2196 } 2197 EXTERN_C_END 2198 2199 #undef __FUNCT__ 2200 #define __FUNCT__ "MatSeqDenseSetLDA" 2201 /*@C 2202 MatSeqDenseSetLDA - Declare the leading dimension of the user-provided array 2203 2204 Input parameter: 2205 + A - the matrix 2206 - lda - the leading dimension 2207 2208 Notes: 2209 This routine is to be used in conjunction with MatSeqDenseSetPreallocation(); 2210 it asserts that the preallocation has a leading dimension (the LDA parameter 2211 of Blas and Lapack fame) larger than M, the first dimension of the matrix. 2212 2213 Level: intermediate 2214 2215 .keywords: dense, matrix, LAPACK, BLAS 2216 2217 .seealso: MatCreate(), MatCreateSeqDense(), MatSeqDenseSetPreallocation(), MatSetMaximumSize() 2218 2219 @*/ 2220 PetscErrorCode MatSeqDenseSetLDA(Mat B,PetscInt lda) 2221 { 2222 Mat_SeqDense *b = (Mat_SeqDense*)B->data; 2223 2224 PetscFunctionBegin; 2225 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); 2226 b->lda = lda; 2227 b->changelda = PETSC_FALSE; 2228 b->Mmax = PetscMax(b->Mmax,lda); 2229 PetscFunctionReturn(0); 2230 } 2231 2232 /*MC 2233 MATSEQDENSE - MATSEQDENSE = "seqdense" - A matrix type to be used for sequential dense matrices. 2234 2235 Options Database Keys: 2236 . -mat_type seqdense - sets the matrix type to "seqdense" during a call to MatSetFromOptions() 2237 2238 Level: beginner 2239 2240 .seealso: MatCreateSeqDense() 2241 2242 M*/ 2243 2244 EXTERN_C_BEGIN 2245 #undef __FUNCT__ 2246 #define __FUNCT__ "MatCreate_SeqDense" 2247 PetscErrorCode MatCreate_SeqDense(Mat B) 2248 { 2249 Mat_SeqDense *b; 2250 PetscErrorCode ierr; 2251 PetscMPIInt size; 2252 2253 PetscFunctionBegin; 2254 ierr = MPI_Comm_size(((PetscObject)B)->comm,&size);CHKERRQ(ierr); 2255 if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Comm must be of size 1"); 2256 2257 ierr = PetscNewLog(B,Mat_SeqDense,&b);CHKERRQ(ierr); 2258 ierr = PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 2259 B->data = (void*)b; 2260 2261 b->pivots = 0; 2262 b->roworiented = PETSC_TRUE; 2263 b->v = 0; 2264 b->changelda = PETSC_FALSE; 2265 2266 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatDenseGetArray_C","MatDenseGetArray_SeqDense",MatDenseGetArray_SeqDense);CHKERRQ(ierr); 2267 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatDenseRestoreArray_C","MatDenseRestoreArray_SeqDense",MatDenseRestoreArray_SeqDense);CHKERRQ(ierr); 2268 2269 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqdense_seqaij_C","MatConvert_SeqDense_SeqAIJ",MatConvert_SeqDense_SeqAIJ);CHKERRQ(ierr); 2270 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_petsc_C", 2271 "MatGetFactor_seqdense_petsc", 2272 MatGetFactor_seqdense_petsc);CHKERRQ(ierr); 2273 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqDenseSetPreallocation_C", 2274 "MatSeqDenseSetPreallocation_SeqDense", 2275 MatSeqDenseSetPreallocation_SeqDense);CHKERRQ(ierr); 2276 2277 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatMatMult_seqaij_seqdense_C", 2278 "MatMatMult_SeqAIJ_SeqDense", 2279 MatMatMult_SeqAIJ_SeqDense);CHKERRQ(ierr); 2280 2281 2282 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatMatMultSymbolic_seqaij_seqdense_C", 2283 "MatMatMultSymbolic_SeqAIJ_SeqDense", 2284 MatMatMultSymbolic_SeqAIJ_SeqDense);CHKERRQ(ierr); 2285 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatMatMultNumeric_seqaij_seqdense_C", 2286 "MatMatMultNumeric_SeqAIJ_SeqDense", 2287 MatMatMultNumeric_SeqAIJ_SeqDense);CHKERRQ(ierr); 2288 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQDENSE);CHKERRQ(ierr); 2289 PetscFunctionReturn(0); 2290 } 2291 EXTERN_C_END 2292