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