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 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix must be factored to solve"); 226 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 227 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 228 ierr = PetscLogFlops(2.0*A->cmap->n*A->cmap->n - A->cmap->n);CHKERRQ(ierr); 229 PetscFunctionReturn(0); 230 } 231 232 #undef __FUNCT__ 233 #define __FUNCT__ "MatMatSolve_SeqDense" 234 PetscErrorCode MatMatSolve_SeqDense(Mat A,Mat B,Mat X) 235 { 236 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 237 PetscErrorCode ierr; 238 PetscScalar *b,*x; 239 PetscInt n; 240 PetscBLASInt nrhs,info,m; 241 PetscBool flg; 242 243 PetscFunctionBegin; 244 ierr = PetscBLASIntCast(A->rmap->n,&m);CHKERRQ(ierr); 245 ierr = PetscObjectTypeCompareAny((PetscObject)B,&flg,MATSEQDENSE,MATMPIDENSE,PETSC_NULL);CHKERRQ(ierr); 246 if (!flg) SETERRQ(((PetscObject)A)->comm,PETSC_ERR_ARG_WRONG,"Matrix B must be MATDENSE matrix"); 247 ierr = PetscObjectTypeCompareAny((PetscObject)X,&flg,MATSEQDENSE,MATMPIDENSE,PETSC_NULL);CHKERRQ(ierr); 248 if (!flg) SETERRQ(((PetscObject)A)->comm,PETSC_ERR_ARG_WRONG,"Matrix X must be MATDENSE matrix"); 249 250 ierr = MatGetSize(B,PETSC_NULL,&n);CHKERRQ(ierr); 251 ierr = PetscBLASIntCast(n,&nrhs);CHKERRQ(ierr); 252 ierr = MatDenseGetArray(B,&b);CHKERRQ(ierr); 253 ierr = MatDenseGetArray(X,&x);CHKERRQ(ierr); 254 255 ierr = PetscMemcpy(x,b,m*nrhs*sizeof(PetscScalar));CHKERRQ(ierr); 256 257 if (A->factortype == MAT_FACTOR_LU) { 258 #if defined(PETSC_MISSING_LAPACK_GETRS) 259 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"GETRS - Lapack routine is unavailable."); 260 #else 261 LAPACKgetrs_("N",&m,&nrhs,mat->v,&mat->lda,mat->pivots,x,&m,&info); 262 if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"GETRS - Bad solve"); 263 #endif 264 } else if (A->factortype == MAT_FACTOR_CHOLESKY) { 265 #if defined(PETSC_MISSING_LAPACK_POTRS) 266 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"POTRS - Lapack routine is unavailable."); 267 #else 268 LAPACKpotrs_("L",&m,&nrhs,mat->v,&mat->lda,x,&m,&info); 269 if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"POTRS Bad solve"); 270 #endif 271 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix must be factored to solve"); 272 273 ierr = MatDenseRestoreArray(B,&b);CHKERRQ(ierr); 274 ierr = MatDenseRestoreArray(X,&x);CHKERRQ(ierr); 275 ierr = PetscLogFlops(nrhs*(2.0*m*m - m));CHKERRQ(ierr); 276 PetscFunctionReturn(0); 277 } 278 279 #undef __FUNCT__ 280 #define __FUNCT__ "MatSolveTranspose_SeqDense" 281 PetscErrorCode MatSolveTranspose_SeqDense(Mat A,Vec xx,Vec yy) 282 { 283 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 284 PetscErrorCode ierr; 285 PetscScalar *x,*y; 286 PetscBLASInt one = 1,info,m; 287 288 PetscFunctionBegin; 289 ierr = PetscBLASIntCast(A->rmap->n,&m);CHKERRQ(ierr); 290 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 291 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 292 ierr = PetscMemcpy(y,x,A->rmap->n*sizeof(PetscScalar));CHKERRQ(ierr); 293 /* assume if pivots exist then use LU; else Cholesky */ 294 if (mat->pivots) { 295 #if defined(PETSC_MISSING_LAPACK_GETRS) 296 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"GETRS - Lapack routine is unavailable."); 297 #else 298 LAPACKgetrs_("T",&m,&one,mat->v,&mat->lda,mat->pivots,y,&m,&info); 299 if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"POTRS - Bad solve"); 300 #endif 301 } else { 302 #if defined(PETSC_MISSING_LAPACK_POTRS) 303 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"POTRS - Lapack routine is unavailable."); 304 #else 305 LAPACKpotrs_("L",&m,&one,mat->v,&mat->lda,y,&m,&info); 306 if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"POTRS - Bad solve"); 307 #endif 308 } 309 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 310 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 311 ierr = PetscLogFlops(2.0*A->cmap->n*A->cmap->n - A->cmap->n);CHKERRQ(ierr); 312 PetscFunctionReturn(0); 313 } 314 315 #undef __FUNCT__ 316 #define __FUNCT__ "MatSolveAdd_SeqDense" 317 PetscErrorCode MatSolveAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy) 318 { 319 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 320 PetscErrorCode ierr; 321 PetscScalar *x,*y,sone = 1.0; 322 Vec tmp = 0; 323 PetscBLASInt one = 1,info,m; 324 325 PetscFunctionBegin; 326 ierr = PetscBLASIntCast(A->rmap->n,&m);CHKERRQ(ierr); 327 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 328 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 329 if (!A->rmap->n || !A->cmap->n) PetscFunctionReturn(0); 330 if (yy == zz) { 331 ierr = VecDuplicate(yy,&tmp);CHKERRQ(ierr); 332 ierr = PetscLogObjectParent(A,tmp);CHKERRQ(ierr); 333 ierr = VecCopy(yy,tmp);CHKERRQ(ierr); 334 } 335 ierr = PetscMemcpy(y,x,A->rmap->n*sizeof(PetscScalar));CHKERRQ(ierr); 336 /* assume if pivots exist then use LU; else Cholesky */ 337 if (mat->pivots) { 338 #if defined(PETSC_MISSING_LAPACK_GETRS) 339 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"GETRS - Lapack routine is unavailable."); 340 #else 341 LAPACKgetrs_("N",&m,&one,mat->v,&mat->lda,mat->pivots,y,&m,&info); 342 if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Bad solve"); 343 #endif 344 } else { 345 #if defined(PETSC_MISSING_LAPACK_POTRS) 346 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"POTRS - Lapack routine is unavailable."); 347 #else 348 LAPACKpotrs_("L",&m,&one,mat->v,&mat->lda,y,&m,&info); 349 if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Bad solve"); 350 #endif 351 } 352 if (tmp) { 353 ierr = VecAXPY(yy,sone,tmp);CHKERRQ(ierr); 354 ierr = VecDestroy(&tmp);CHKERRQ(ierr); 355 } else { 356 ierr = VecAXPY(yy,sone,zz);CHKERRQ(ierr); 357 } 358 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 359 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 360 ierr = PetscLogFlops(2.0*A->cmap->n*A->cmap->n);CHKERRQ(ierr); 361 PetscFunctionReturn(0); 362 } 363 364 #undef __FUNCT__ 365 #define __FUNCT__ "MatSolveTransposeAdd_SeqDense" 366 PetscErrorCode MatSolveTransposeAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy) 367 { 368 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 369 PetscErrorCode ierr; 370 PetscScalar *x,*y,sone = 1.0; 371 Vec tmp; 372 PetscBLASInt one = 1,info,m; 373 374 PetscFunctionBegin; 375 ierr = PetscBLASIntCast(A->rmap->n,&m);CHKERRQ(ierr); 376 if (!A->rmap->n || !A->cmap->n) PetscFunctionReturn(0); 377 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 378 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 379 if (yy == zz) { 380 ierr = VecDuplicate(yy,&tmp);CHKERRQ(ierr); 381 ierr = PetscLogObjectParent(A,tmp);CHKERRQ(ierr); 382 ierr = VecCopy(yy,tmp);CHKERRQ(ierr); 383 } 384 ierr = PetscMemcpy(y,x,A->rmap->n*sizeof(PetscScalar));CHKERRQ(ierr); 385 /* assume if pivots exist then use LU; else Cholesky */ 386 if (mat->pivots) { 387 #if defined(PETSC_MISSING_LAPACK_GETRS) 388 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"GETRS - Lapack routine is unavailable."); 389 #else 390 LAPACKgetrs_("T",&m,&one,mat->v,&mat->lda,mat->pivots,y,&m,&info); 391 if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Bad solve"); 392 #endif 393 } else { 394 #if defined(PETSC_MISSING_LAPACK_POTRS) 395 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"POTRS - Lapack routine is unavailable."); 396 #else 397 LAPACKpotrs_("L",&m,&one,mat->v,&mat->lda,y,&m,&info); 398 if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Bad solve"); 399 #endif 400 } 401 if (tmp) { 402 ierr = VecAXPY(yy,sone,tmp);CHKERRQ(ierr); 403 ierr = VecDestroy(&tmp);CHKERRQ(ierr); 404 } else { 405 ierr = VecAXPY(yy,sone,zz);CHKERRQ(ierr); 406 } 407 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 408 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 409 ierr = PetscLogFlops(2.0*A->cmap->n*A->cmap->n);CHKERRQ(ierr); 410 PetscFunctionReturn(0); 411 } 412 413 /* ---------------------------------------------------------------*/ 414 /* COMMENT: I have chosen to hide row permutation in the pivots, 415 rather than put it in the Mat->row slot.*/ 416 #undef __FUNCT__ 417 #define __FUNCT__ "MatLUFactor_SeqDense" 418 PetscErrorCode MatLUFactor_SeqDense(Mat A,IS row,IS col,const MatFactorInfo *minfo) 419 { 420 #if defined(PETSC_MISSING_LAPACK_GETRF) 421 PetscFunctionBegin; 422 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"GETRF - Lapack routine is unavailable."); 423 #else 424 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 425 PetscErrorCode ierr; 426 PetscBLASInt n,m,info; 427 428 PetscFunctionBegin; 429 ierr = PetscBLASIntCast(A->cmap->n,&n);CHKERRQ(ierr); 430 ierr = PetscBLASIntCast(A->rmap->n,&m);CHKERRQ(ierr); 431 if (!mat->pivots) { 432 ierr = PetscMalloc((A->rmap->n+1)*sizeof(PetscBLASInt),&mat->pivots);CHKERRQ(ierr); 433 ierr = PetscLogObjectMemory(A,A->rmap->n*sizeof(PetscBLASInt));CHKERRQ(ierr); 434 } 435 if (!A->rmap->n || !A->cmap->n) PetscFunctionReturn(0); 436 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr); 437 LAPACKgetrf_(&m,&n,mat->v,&mat->lda,mat->pivots,&info); 438 ierr = PetscFPTrapPop();CHKERRQ(ierr); 439 440 if (info<0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Bad argument to LU factorization"); 441 if (info>0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_MAT_LU_ZRPVT,"Bad LU factorization"); 442 A->ops->solve = MatSolve_SeqDense; 443 A->ops->solvetranspose = MatSolveTranspose_SeqDense; 444 A->ops->solveadd = MatSolveAdd_SeqDense; 445 A->ops->solvetransposeadd = MatSolveTransposeAdd_SeqDense; 446 A->factortype = MAT_FACTOR_LU; 447 448 ierr = PetscLogFlops((2.0*A->cmap->n*A->cmap->n*A->cmap->n)/3);CHKERRQ(ierr); 449 #endif 450 PetscFunctionReturn(0); 451 } 452 453 #undef __FUNCT__ 454 #define __FUNCT__ "MatCholeskyFactor_SeqDense" 455 PetscErrorCode MatCholeskyFactor_SeqDense(Mat A,IS perm,const MatFactorInfo *factinfo) 456 { 457 #if defined(PETSC_MISSING_LAPACK_POTRF) 458 PetscFunctionBegin; 459 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"POTRF - Lapack routine is unavailable."); 460 #else 461 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 462 PetscErrorCode ierr; 463 PetscBLASInt info,n; 464 465 PetscFunctionBegin; 466 ierr = PetscBLASIntCast(A->cmap->n,&n);CHKERRQ(ierr); 467 ierr = PetscFree(mat->pivots);CHKERRQ(ierr); 468 469 if (!A->rmap->n || !A->cmap->n) PetscFunctionReturn(0); 470 LAPACKpotrf_("L",&n,mat->v,&mat->lda,&info); 471 if (info) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_MAT_CH_ZRPVT,"Bad factorization: zero pivot in row %D",(PetscInt)info-1); 472 A->ops->solve = MatSolve_SeqDense; 473 A->ops->solvetranspose = MatSolveTranspose_SeqDense; 474 A->ops->solveadd = MatSolveAdd_SeqDense; 475 A->ops->solvetransposeadd = MatSolveTransposeAdd_SeqDense; 476 A->factortype = MAT_FACTOR_CHOLESKY; 477 478 ierr = PetscLogFlops((A->cmap->n*A->cmap->n*A->cmap->n)/3.0);CHKERRQ(ierr); 479 #endif 480 PetscFunctionReturn(0); 481 } 482 483 484 #undef __FUNCT__ 485 #define __FUNCT__ "MatCholeskyFactorNumeric_SeqDense" 486 PetscErrorCode MatCholeskyFactorNumeric_SeqDense(Mat fact,Mat A,const MatFactorInfo *info_dummy) 487 { 488 PetscErrorCode ierr; 489 MatFactorInfo info; 490 491 PetscFunctionBegin; 492 info.fill = 1.0; 493 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 986 col_lens[0] = MAT_FILE_CLASSID; 987 col_lens[1] = m; 988 col_lens[2] = n; 989 col_lens[3] = MATRIX_BINARY_FORMAT_DENSE; 990 991 ierr = PetscBinaryWrite(fd,col_lens,4,PETSC_INT,PETSC_TRUE);CHKERRQ(ierr); 992 ierr = PetscFree(col_lens);CHKERRQ(ierr); 993 994 /* write out matrix, by rows */ 995 ierr = PetscMalloc((m*n+1)*sizeof(PetscScalar),&vals);CHKERRQ(ierr); 996 v = a->v; 997 for (j=0; j<n; j++) { 998 for (i=0; i<m; i++) { 999 vals[j + i*n] = *v++; 1000 } 1001 } 1002 ierr = PetscBinaryWrite(fd,vals,n*m,PETSC_SCALAR,PETSC_FALSE);CHKERRQ(ierr); 1003 ierr = PetscFree(vals);CHKERRQ(ierr); 1004 } else { 1005 ierr = PetscMalloc((4+nz)*sizeof(PetscInt),&col_lens);CHKERRQ(ierr); 1006 1007 col_lens[0] = MAT_FILE_CLASSID; 1008 col_lens[1] = m; 1009 col_lens[2] = n; 1010 col_lens[3] = nz; 1011 1012 /* store lengths of each row and write (including header) to file */ 1013 for (i=0; i<m; i++) col_lens[4+i] = n; 1014 ierr = PetscBinaryWrite(fd,col_lens,4+m,PETSC_INT,PETSC_TRUE);CHKERRQ(ierr); 1015 1016 /* Possibly should write in smaller increments, not whole matrix at once? */ 1017 /* store column indices (zero start index) */ 1018 ict = 0; 1019 for (i=0; i<m; i++) { 1020 for (j=0; j<n; j++) col_lens[ict++] = j; 1021 } 1022 ierr = PetscBinaryWrite(fd,col_lens,nz,PETSC_INT,PETSC_FALSE);CHKERRQ(ierr); 1023 ierr = PetscFree(col_lens);CHKERRQ(ierr); 1024 1025 /* store nonzero values */ 1026 ierr = PetscMalloc((nz+1)*sizeof(PetscScalar),&anonz);CHKERRQ(ierr); 1027 ict = 0; 1028 for (i=0; i<m; i++) { 1029 v = a->v + i; 1030 for (j=0; j<n; j++) { 1031 anonz[ict++] = *v; v += a->lda; 1032 } 1033 } 1034 ierr = PetscBinaryWrite(fd,anonz,nz,PETSC_SCALAR,PETSC_FALSE);CHKERRQ(ierr); 1035 ierr = PetscFree(anonz);CHKERRQ(ierr); 1036 } 1037 PetscFunctionReturn(0); 1038 } 1039 1040 #undef __FUNCT__ 1041 #define __FUNCT__ "MatView_SeqDense_Draw_Zoom" 1042 PetscErrorCode MatView_SeqDense_Draw_Zoom(PetscDraw draw,void *Aa) 1043 { 1044 Mat A = (Mat) Aa; 1045 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1046 PetscErrorCode ierr; 1047 PetscInt m = A->rmap->n,n = A->cmap->n,color,i,j; 1048 PetscScalar *v = a->v; 1049 PetscViewer viewer; 1050 PetscDraw popup; 1051 PetscReal xl,yl,xr,yr,x_l,x_r,y_l,y_r,scale,maxv = 0.0; 1052 PetscViewerFormat format; 1053 1054 PetscFunctionBegin; 1055 ierr = PetscObjectQuery((PetscObject)A,"Zoomviewer",(PetscObject*)&viewer);CHKERRQ(ierr); 1056 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 1057 ierr = PetscDrawGetCoordinates(draw,&xl,&yl,&xr,&yr);CHKERRQ(ierr); 1058 1059 /* Loop over matrix elements drawing boxes */ 1060 if (format != PETSC_VIEWER_DRAW_CONTOUR) { 1061 /* Blue for negative and Red for positive */ 1062 color = PETSC_DRAW_BLUE; 1063 for (j = 0; j < n; j++) { 1064 x_l = j; 1065 x_r = x_l + 1.0; 1066 for (i = 0; i < m; i++) { 1067 y_l = m - i - 1.0; 1068 y_r = y_l + 1.0; 1069 if (PetscRealPart(v[j*m+i]) > 0.) { 1070 color = PETSC_DRAW_RED; 1071 } else if (PetscRealPart(v[j*m+i]) < 0.) { 1072 color = PETSC_DRAW_BLUE; 1073 } else { 1074 continue; 1075 } 1076 ierr = PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);CHKERRQ(ierr); 1077 } 1078 } 1079 } else { 1080 /* use contour shading to indicate magnitude of values */ 1081 /* first determine max of all nonzero values */ 1082 for (i = 0; i < m*n; i++) { 1083 if (PetscAbsScalar(v[i]) > maxv) maxv = PetscAbsScalar(v[i]); 1084 } 1085 scale = (245.0 - PETSC_DRAW_BASIC_COLORS)/maxv; 1086 ierr = PetscDrawGetPopup(draw,&popup);CHKERRQ(ierr); 1087 if (popup) {ierr = PetscDrawScalePopup(popup,0.0,maxv);CHKERRQ(ierr);} 1088 for (j = 0; j < n; j++) { 1089 x_l = j; 1090 x_r = x_l + 1.0; 1091 for (i = 0; i < m; i++) { 1092 y_l = m - i - 1.0; 1093 y_r = y_l + 1.0; 1094 color = PETSC_DRAW_BASIC_COLORS + (int)(scale*PetscAbsScalar(v[j*m+i])); 1095 ierr = PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);CHKERRQ(ierr); 1096 } 1097 } 1098 } 1099 PetscFunctionReturn(0); 1100 } 1101 1102 #undef __FUNCT__ 1103 #define __FUNCT__ "MatView_SeqDense_Draw" 1104 PetscErrorCode MatView_SeqDense_Draw(Mat A,PetscViewer viewer) 1105 { 1106 PetscDraw draw; 1107 PetscBool isnull; 1108 PetscReal xr,yr,xl,yl,h,w; 1109 PetscErrorCode ierr; 1110 1111 PetscFunctionBegin; 1112 ierr = PetscViewerDrawGetDraw(viewer,0,&draw);CHKERRQ(ierr); 1113 ierr = PetscDrawIsNull(draw,&isnull);CHKERRQ(ierr); 1114 if (isnull) PetscFunctionReturn(0); 1115 1116 ierr = PetscObjectCompose((PetscObject)A,"Zoomviewer",(PetscObject)viewer);CHKERRQ(ierr); 1117 xr = A->cmap->n; yr = A->rmap->n; h = yr/10.0; w = xr/10.0; 1118 xr += w; yr += h; xl = -w; yl = -h; 1119 ierr = PetscDrawSetCoordinates(draw,xl,yl,xr,yr);CHKERRQ(ierr); 1120 ierr = PetscDrawZoom(draw,MatView_SeqDense_Draw_Zoom,A);CHKERRQ(ierr); 1121 ierr = PetscObjectCompose((PetscObject)A,"Zoomviewer",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 } 1144 PetscFunctionReturn(0); 1145 } 1146 1147 #undef __FUNCT__ 1148 #define __FUNCT__ "MatDestroy_SeqDense" 1149 PetscErrorCode MatDestroy_SeqDense(Mat mat) 1150 { 1151 Mat_SeqDense *l = (Mat_SeqDense*)mat->data; 1152 PetscErrorCode ierr; 1153 1154 PetscFunctionBegin; 1155 #if defined(PETSC_USE_LOG) 1156 PetscLogObjectState((PetscObject)mat,"Rows %D Cols %D",mat->rmap->n,mat->cmap->n); 1157 #endif 1158 ierr = PetscFree(l->pivots);CHKERRQ(ierr); 1159 if (!l->user_alloc) {ierr = PetscFree(l->v);CHKERRQ(ierr);} 1160 ierr = PetscFree(mat->data);CHKERRQ(ierr); 1161 1162 ierr = PetscObjectChangeTypeName((PetscObject)mat,0);CHKERRQ(ierr); 1163 ierr = PetscObjectComposeFunctionDynamic((PetscObject)mat,"MatDenseGetArray_C","",PETSC_NULL);CHKERRQ(ierr); 1164 ierr = PetscObjectComposeFunctionDynamic((PetscObject)mat,"MatDenseRestoreArray_C","",PETSC_NULL);CHKERRQ(ierr); 1165 ierr = PetscObjectComposeFunctionDynamic((PetscObject)mat,"MatSeqDenseSetPreallocation_C","",PETSC_NULL);CHKERRQ(ierr); 1166 ierr = PetscObjectComposeFunctionDynamic((PetscObject)mat,"MatMatMult_seqaij_seqdense_C","",PETSC_NULL);CHKERRQ(ierr); 1167 ierr = PetscObjectComposeFunctionDynamic((PetscObject)mat,"MatMatMultSymbolic_seqaij_seqdense_C","",PETSC_NULL);CHKERRQ(ierr); 1168 ierr = PetscObjectComposeFunctionDynamic((PetscObject)mat,"MatMatMultNumeric_seqaij_seqdense_C","",PETSC_NULL);CHKERRQ(ierr); 1169 PetscFunctionReturn(0); 1170 } 1171 1172 #undef __FUNCT__ 1173 #define __FUNCT__ "MatTranspose_SeqDense" 1174 PetscErrorCode MatTranspose_SeqDense(Mat A,MatReuse reuse,Mat *matout) 1175 { 1176 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1177 PetscErrorCode ierr; 1178 PetscInt k,j,m,n,M; 1179 PetscScalar *v,tmp; 1180 1181 PetscFunctionBegin; 1182 v = mat->v; m = A->rmap->n; M = mat->lda; n = A->cmap->n; 1183 if (reuse == MAT_REUSE_MATRIX && *matout == A) { /* in place transpose */ 1184 if (m != n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Can not transpose non-square matrix in place"); 1185 else { 1186 for (j=0; j<m; j++) { 1187 for (k=0; k<j; k++) { 1188 tmp = v[j + k*M]; 1189 v[j + k*M] = v[k + j*M]; 1190 v[k + j*M] = tmp; 1191 } 1192 } 1193 } 1194 } else { /* out-of-place transpose */ 1195 Mat tmat; 1196 Mat_SeqDense *tmatd; 1197 PetscScalar *v2; 1198 1199 if (reuse == MAT_INITIAL_MATRIX) { 1200 ierr = MatCreate(((PetscObject)A)->comm,&tmat);CHKERRQ(ierr); 1201 ierr = MatSetSizes(tmat,A->cmap->n,A->rmap->n,A->cmap->n,A->rmap->n);CHKERRQ(ierr); 1202 ierr = MatSetType(tmat,((PetscObject)A)->type_name);CHKERRQ(ierr); 1203 ierr = MatSeqDenseSetPreallocation(tmat,PETSC_NULL);CHKERRQ(ierr); 1204 } else { 1205 tmat = *matout; 1206 } 1207 tmatd = (Mat_SeqDense*)tmat->data; 1208 v = mat->v; v2 = tmatd->v; 1209 for (j=0; j<n; j++) { 1210 for (k=0; k<m; k++) v2[j + k*n] = v[k + j*M]; 1211 } 1212 ierr = MatAssemblyBegin(tmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1213 ierr = MatAssemblyEnd(tmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1214 *matout = tmat; 1215 } 1216 PetscFunctionReturn(0); 1217 } 1218 1219 #undef __FUNCT__ 1220 #define __FUNCT__ "MatEqual_SeqDense" 1221 PetscErrorCode MatEqual_SeqDense(Mat A1,Mat A2,PetscBool *flg) 1222 { 1223 Mat_SeqDense *mat1 = (Mat_SeqDense*)A1->data; 1224 Mat_SeqDense *mat2 = (Mat_SeqDense*)A2->data; 1225 PetscInt i,j; 1226 PetscScalar *v1,*v2; 1227 1228 PetscFunctionBegin; 1229 if (A1->rmap->n != A2->rmap->n) {*flg = PETSC_FALSE; PetscFunctionReturn(0);} 1230 if (A1->cmap->n != A2->cmap->n) {*flg = PETSC_FALSE; PetscFunctionReturn(0);} 1231 for (i=0; i<A1->rmap->n; i++) { 1232 v1 = mat1->v+i; v2 = mat2->v+i; 1233 for (j=0; j<A1->cmap->n; j++) { 1234 if (*v1 != *v2) {*flg = PETSC_FALSE; PetscFunctionReturn(0);} 1235 v1 += mat1->lda; v2 += mat2->lda; 1236 } 1237 } 1238 *flg = PETSC_TRUE; 1239 PetscFunctionReturn(0); 1240 } 1241 1242 #undef __FUNCT__ 1243 #define __FUNCT__ "MatGetDiagonal_SeqDense" 1244 PetscErrorCode MatGetDiagonal_SeqDense(Mat A,Vec v) 1245 { 1246 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1247 PetscErrorCode ierr; 1248 PetscInt i,n,len; 1249 PetscScalar *x,zero = 0.0; 1250 1251 PetscFunctionBegin; 1252 ierr = VecSet(v,zero);CHKERRQ(ierr); 1253 ierr = VecGetSize(v,&n);CHKERRQ(ierr); 1254 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 1255 len = PetscMin(A->rmap->n,A->cmap->n); 1256 if (n != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming mat and vec"); 1257 for (i=0; i<len; i++) { 1258 x[i] = mat->v[i*mat->lda + i]; 1259 } 1260 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 1261 PetscFunctionReturn(0); 1262 } 1263 1264 #undef __FUNCT__ 1265 #define __FUNCT__ "MatDiagonalScale_SeqDense" 1266 PetscErrorCode MatDiagonalScale_SeqDense(Mat A,Vec ll,Vec rr) 1267 { 1268 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1269 PetscScalar *l,*r,x,*v; 1270 PetscErrorCode ierr; 1271 PetscInt i,j,m = A->rmap->n,n = A->cmap->n; 1272 1273 PetscFunctionBegin; 1274 if (ll) { 1275 ierr = VecGetSize(ll,&m);CHKERRQ(ierr); 1276 ierr = VecGetArray(ll,&l);CHKERRQ(ierr); 1277 if (m != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Left scaling vec wrong size"); 1278 for (i=0; i<m; i++) { 1279 x = l[i]; 1280 v = mat->v + i; 1281 for (j=0; j<n; j++) { (*v) *= x; v+= m;} 1282 } 1283 ierr = VecRestoreArray(ll,&l);CHKERRQ(ierr); 1284 ierr = PetscLogFlops(n*m);CHKERRQ(ierr); 1285 } 1286 if (rr) { 1287 ierr = VecGetSize(rr,&n);CHKERRQ(ierr); 1288 ierr = VecGetArray(rr,&r);CHKERRQ(ierr); 1289 if (n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Right scaling vec wrong size"); 1290 for (i=0; i<n; i++) { 1291 x = r[i]; 1292 v = mat->v + i*m; 1293 for (j=0; j<m; j++) (*v++) *= x; 1294 } 1295 ierr = VecRestoreArray(rr,&r);CHKERRQ(ierr); 1296 ierr = PetscLogFlops(n*m);CHKERRQ(ierr); 1297 } 1298 PetscFunctionReturn(0); 1299 } 1300 1301 #undef __FUNCT__ 1302 #define __FUNCT__ "MatNorm_SeqDense" 1303 PetscErrorCode MatNorm_SeqDense(Mat A,NormType type,PetscReal *nrm) 1304 { 1305 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1306 PetscScalar *v = mat->v; 1307 PetscReal sum = 0.0; 1308 PetscInt lda =mat->lda,m=A->rmap->n,i,j; 1309 PetscErrorCode ierr; 1310 1311 PetscFunctionBegin; 1312 if (type == NORM_FROBENIUS) { 1313 if (lda>m) { 1314 for (j=0; j<A->cmap->n; j++) { 1315 v = mat->v+j*lda; 1316 for (i=0; i<m; i++) { 1317 sum += PetscRealPart(PetscConj(*v)*(*v)); v++; 1318 } 1319 } 1320 } else { 1321 for (i=0; i<A->cmap->n*A->rmap->n; i++) { 1322 sum += PetscRealPart(PetscConj(*v)*(*v)); v++; 1323 } 1324 } 1325 *nrm = PetscSqrtReal(sum); 1326 ierr = PetscLogFlops(2.0*A->cmap->n*A->rmap->n);CHKERRQ(ierr); 1327 } else if (type == NORM_1) { 1328 *nrm = 0.0; 1329 for (j=0; j<A->cmap->n; j++) { 1330 v = mat->v + j*mat->lda; 1331 sum = 0.0; 1332 for (i=0; i<A->rmap->n; i++) { 1333 sum += PetscAbsScalar(*v); v++; 1334 } 1335 if (sum > *nrm) *nrm = sum; 1336 } 1337 ierr = PetscLogFlops(A->cmap->n*A->rmap->n);CHKERRQ(ierr); 1338 } else if (type == NORM_INFINITY) { 1339 *nrm = 0.0; 1340 for (j=0; j<A->rmap->n; j++) { 1341 v = mat->v + j; 1342 sum = 0.0; 1343 for (i=0; i<A->cmap->n; i++) { 1344 sum += PetscAbsScalar(*v); v += mat->lda; 1345 } 1346 if (sum > *nrm) *nrm = sum; 1347 } 1348 ierr = PetscLogFlops(A->cmap->n*A->rmap->n);CHKERRQ(ierr); 1349 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No two norm"); 1350 PetscFunctionReturn(0); 1351 } 1352 1353 #undef __FUNCT__ 1354 #define __FUNCT__ "MatSetOption_SeqDense" 1355 PetscErrorCode MatSetOption_SeqDense(Mat A,MatOption op,PetscBool flg) 1356 { 1357 Mat_SeqDense *aij = (Mat_SeqDense*)A->data; 1358 PetscErrorCode ierr; 1359 1360 PetscFunctionBegin; 1361 switch (op) { 1362 case MAT_ROW_ORIENTED: 1363 aij->roworiented = flg; 1364 break; 1365 case MAT_NEW_NONZERO_LOCATIONS: 1366 case MAT_NEW_NONZERO_LOCATION_ERR: 1367 case MAT_NEW_NONZERO_ALLOCATION_ERR: 1368 case MAT_NEW_DIAGONALS: 1369 case MAT_KEEP_NONZERO_PATTERN: 1370 case MAT_IGNORE_OFF_PROC_ENTRIES: 1371 case MAT_USE_HASH_TABLE: 1372 case MAT_IGNORE_LOWER_TRIANGULAR: 1373 ierr = PetscInfo1(A,"Option %s ignored\n",MatOptions[op]);CHKERRQ(ierr); 1374 break; 1375 case MAT_SPD: 1376 case MAT_SYMMETRIC: 1377 case MAT_STRUCTURALLY_SYMMETRIC: 1378 case MAT_HERMITIAN: 1379 case MAT_SYMMETRY_ETERNAL: 1380 /* These options are handled directly by MatSetOption() */ 1381 break; 1382 default: 1383 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"unknown option %s",MatOptions[op]); 1384 } 1385 PetscFunctionReturn(0); 1386 } 1387 1388 #undef __FUNCT__ 1389 #define __FUNCT__ "MatZeroEntries_SeqDense" 1390 PetscErrorCode MatZeroEntries_SeqDense(Mat A) 1391 { 1392 Mat_SeqDense *l = (Mat_SeqDense*)A->data; 1393 PetscErrorCode ierr; 1394 PetscInt lda=l->lda,m=A->rmap->n,j; 1395 1396 PetscFunctionBegin; 1397 if (lda>m) { 1398 for (j=0; j<A->cmap->n; j++) { 1399 ierr = PetscMemzero(l->v+j*lda,m*sizeof(PetscScalar));CHKERRQ(ierr); 1400 } 1401 } else { 1402 ierr = PetscMemzero(l->v,A->rmap->n*A->cmap->n*sizeof(PetscScalar));CHKERRQ(ierr); 1403 } 1404 PetscFunctionReturn(0); 1405 } 1406 1407 #undef __FUNCT__ 1408 #define __FUNCT__ "MatZeroRows_SeqDense" 1409 PetscErrorCode MatZeroRows_SeqDense(Mat A,PetscInt N,const PetscInt rows[],PetscScalar diag,Vec x,Vec b) 1410 { 1411 PetscErrorCode ierr; 1412 Mat_SeqDense *l = (Mat_SeqDense*)A->data; 1413 PetscInt m = l->lda, n = A->cmap->n, i,j; 1414 PetscScalar *slot,*bb; 1415 const PetscScalar *xx; 1416 1417 PetscFunctionBegin; 1418 #if defined(PETSC_USE_DEBUG) 1419 for (i=0; i<N; i++) { 1420 if (rows[i] < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative row requested to be zeroed"); 1421 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); 1422 } 1423 #endif 1424 1425 /* fix right hand side if needed */ 1426 if (x && b) { 1427 ierr = VecGetArrayRead(x,&xx);CHKERRQ(ierr); 1428 ierr = VecGetArray(b,&bb);CHKERRQ(ierr); 1429 for (i=0; i<N; i++) bb[rows[i]] = diag*xx[rows[i]]; 1430 ierr = VecRestoreArrayRead(x,&xx);CHKERRQ(ierr); 1431 ierr = VecRestoreArray(b,&bb);CHKERRQ(ierr); 1432 } 1433 1434 for (i=0; i<N; i++) { 1435 slot = l->v + rows[i]; 1436 for (j=0; j<n; j++) { *slot = 0.0; slot += m;} 1437 } 1438 if (diag != 0.0) { 1439 if (A->rmap->n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only coded for square matrices"); 1440 for (i=0; i<N; i++) { 1441 slot = l->v + (m+1)*rows[i]; 1442 *slot = diag; 1443 } 1444 } 1445 PetscFunctionReturn(0); 1446 } 1447 1448 #undef __FUNCT__ 1449 #define __FUNCT__ "MatDenseGetArray_SeqDense" 1450 PetscErrorCode MatDenseGetArray_SeqDense(Mat A,PetscScalar *array[]) 1451 { 1452 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1453 1454 PetscFunctionBegin; 1455 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"); 1456 *array = mat->v; 1457 PetscFunctionReturn(0); 1458 } 1459 1460 #undef __FUNCT__ 1461 #define __FUNCT__ "MatDenseRestoreArray_SeqDense" 1462 PetscErrorCode MatDenseRestoreArray_SeqDense(Mat A,PetscScalar *array[]) 1463 { 1464 PetscFunctionBegin; 1465 *array = 0; /* user cannot accidently use the array later */ 1466 PetscFunctionReturn(0); 1467 } 1468 1469 #undef __FUNCT__ 1470 #define __FUNCT__ "MatDenseGetArray" 1471 /*@C 1472 MatDenseGetArray - gives access to the array where the data for a SeqDense matrix is stored 1473 1474 Not Collective 1475 1476 Input Parameter: 1477 . mat - a MATSEQDENSE matrix 1478 1479 Output Parameter: 1480 . array - pointer to the data 1481 1482 Level: intermediate 1483 1484 .seealso: MatDenseRestoreArray() 1485 @*/ 1486 PetscErrorCode MatDenseGetArray(Mat A,PetscScalar **array) 1487 { 1488 PetscErrorCode ierr; 1489 1490 PetscFunctionBegin; 1491 ierr = PetscUseMethod(A,"MatDenseGetArray_C",(Mat,PetscScalar**),(A,array));CHKERRQ(ierr); 1492 PetscFunctionReturn(0); 1493 } 1494 1495 #undef __FUNCT__ 1496 #define __FUNCT__ "MatDenseRestoreArray" 1497 /*@C 1498 MatDenseRestoreArray - returns access to the array where the data for a SeqDense matrix is stored obtained by MatDenseGetArray() 1499 1500 Not Collective 1501 1502 Input Parameters: 1503 . mat - a MATSEQDENSE matrix 1504 . array - pointer to the data 1505 1506 Level: intermediate 1507 1508 .seealso: MatDenseGetArray() 1509 @*/ 1510 PetscErrorCode MatDenseRestoreArray(Mat A,PetscScalar **array) 1511 { 1512 PetscErrorCode ierr; 1513 1514 PetscFunctionBegin; 1515 ierr = PetscUseMethod(A,"MatDenseRestoreArray_C",(Mat,PetscScalar**),(A,array));CHKERRQ(ierr); 1516 PetscFunctionReturn(0); 1517 } 1518 1519 #undef __FUNCT__ 1520 #define __FUNCT__ "MatGetSubMatrix_SeqDense" 1521 static PetscErrorCode MatGetSubMatrix_SeqDense(Mat A,IS isrow,IS iscol,PetscInt cs,MatReuse scall,Mat *B) 1522 { 1523 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1524 PetscErrorCode ierr; 1525 PetscInt i,j,nrows,ncols; 1526 const PetscInt *irow,*icol; 1527 PetscScalar *av,*bv,*v = mat->v; 1528 Mat newmat; 1529 1530 PetscFunctionBegin; 1531 ierr = ISGetIndices(isrow,&irow);CHKERRQ(ierr); 1532 ierr = ISGetIndices(iscol,&icol);CHKERRQ(ierr); 1533 ierr = ISGetLocalSize(isrow,&nrows);CHKERRQ(ierr); 1534 ierr = ISGetLocalSize(iscol,&ncols);CHKERRQ(ierr); 1535 1536 /* Check submatrixcall */ 1537 if (scall == MAT_REUSE_MATRIX) { 1538 PetscInt n_cols,n_rows; 1539 ierr = MatGetSize(*B,&n_rows,&n_cols);CHKERRQ(ierr); 1540 if (n_rows != nrows || n_cols != ncols) { 1541 /* resize the result matrix to match number of requested rows/columns */ 1542 ierr = MatSetSizes(*B,nrows,ncols,nrows,ncols);CHKERRQ(ierr); 1543 } 1544 newmat = *B; 1545 } else { 1546 /* Create and fill new matrix */ 1547 ierr = MatCreate(((PetscObject)A)->comm,&newmat);CHKERRQ(ierr); 1548 ierr = MatSetSizes(newmat,nrows,ncols,nrows,ncols);CHKERRQ(ierr); 1549 ierr = MatSetType(newmat,((PetscObject)A)->type_name);CHKERRQ(ierr); 1550 ierr = MatSeqDenseSetPreallocation(newmat,PETSC_NULL);CHKERRQ(ierr); 1551 } 1552 1553 /* Now extract the data pointers and do the copy,column at a time */ 1554 bv = ((Mat_SeqDense*)newmat->data)->v; 1555 1556 for (i=0; i<ncols; i++) { 1557 av = v + mat->lda*icol[i]; 1558 for (j=0; j<nrows; j++) *bv++ = av[irow[j]]; 1559 } 1560 1561 /* Assemble the matrices so that the correct flags are set */ 1562 ierr = MatAssemblyBegin(newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1563 ierr = MatAssemblyEnd(newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1564 1565 /* Free work space */ 1566 ierr = ISRestoreIndices(isrow,&irow);CHKERRQ(ierr); 1567 ierr = ISRestoreIndices(iscol,&icol);CHKERRQ(ierr); 1568 *B = newmat; 1569 PetscFunctionReturn(0); 1570 } 1571 1572 #undef __FUNCT__ 1573 #define __FUNCT__ "MatGetSubMatrices_SeqDense" 1574 PetscErrorCode MatGetSubMatrices_SeqDense(Mat A,PetscInt n,const IS irow[],const IS icol[],MatReuse scall,Mat *B[]) 1575 { 1576 PetscErrorCode ierr; 1577 PetscInt i; 1578 1579 PetscFunctionBegin; 1580 if (scall == MAT_INITIAL_MATRIX) { 1581 ierr = PetscMalloc((n+1)*sizeof(Mat),B);CHKERRQ(ierr); 1582 } 1583 1584 for (i=0; i<n; i++) { 1585 ierr = MatGetSubMatrix_SeqDense(A,irow[i],icol[i],PETSC_DECIDE,scall,&(*B)[i]);CHKERRQ(ierr); 1586 } 1587 PetscFunctionReturn(0); 1588 } 1589 1590 #undef __FUNCT__ 1591 #define __FUNCT__ "MatAssemblyBegin_SeqDense" 1592 PetscErrorCode MatAssemblyBegin_SeqDense(Mat mat,MatAssemblyType mode) 1593 { 1594 PetscFunctionBegin; 1595 PetscFunctionReturn(0); 1596 } 1597 1598 #undef __FUNCT__ 1599 #define __FUNCT__ "MatAssemblyEnd_SeqDense" 1600 PetscErrorCode MatAssemblyEnd_SeqDense(Mat mat,MatAssemblyType mode) 1601 { 1602 PetscFunctionBegin; 1603 PetscFunctionReturn(0); 1604 } 1605 1606 #undef __FUNCT__ 1607 #define __FUNCT__ "MatCopy_SeqDense" 1608 PetscErrorCode MatCopy_SeqDense(Mat A,Mat B,MatStructure str) 1609 { 1610 Mat_SeqDense *a = (Mat_SeqDense*)A->data,*b = (Mat_SeqDense*)B->data; 1611 PetscErrorCode ierr; 1612 PetscInt lda1=a->lda,lda2=b->lda, m=A->rmap->n,n=A->cmap->n, j; 1613 1614 PetscFunctionBegin; 1615 /* If the two matrices don't have the same copy implementation, they aren't compatible for fast copy. */ 1616 if (A->ops->copy != B->ops->copy) { 1617 ierr = MatCopy_Basic(A,B,str);CHKERRQ(ierr); 1618 PetscFunctionReturn(0); 1619 } 1620 if (m != B->rmap->n || n != B->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"size(B) != size(A)"); 1621 if (lda1>m || lda2>m) { 1622 for (j=0; j<n; j++) { 1623 ierr = PetscMemcpy(b->v+j*lda2,a->v+j*lda1,m*sizeof(PetscScalar));CHKERRQ(ierr); 1624 } 1625 } else { 1626 ierr = PetscMemcpy(b->v,a->v,A->rmap->n*A->cmap->n*sizeof(PetscScalar));CHKERRQ(ierr); 1627 } 1628 PetscFunctionReturn(0); 1629 } 1630 1631 #undef __FUNCT__ 1632 #define __FUNCT__ "MatSetUp_SeqDense" 1633 PetscErrorCode MatSetUp_SeqDense(Mat A) 1634 { 1635 PetscErrorCode ierr; 1636 1637 PetscFunctionBegin; 1638 ierr = MatSeqDenseSetPreallocation(A,0);CHKERRQ(ierr); 1639 PetscFunctionReturn(0); 1640 } 1641 1642 #undef __FUNCT__ 1643 #define __FUNCT__ "MatConjugate_SeqDense" 1644 static PetscErrorCode MatConjugate_SeqDense(Mat A) 1645 { 1646 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1647 PetscInt i,nz = A->rmap->n*A->cmap->n; 1648 PetscScalar *aa = a->v; 1649 1650 PetscFunctionBegin; 1651 for (i=0; i<nz; i++) aa[i] = PetscConj(aa[i]); 1652 PetscFunctionReturn(0); 1653 } 1654 1655 #undef __FUNCT__ 1656 #define __FUNCT__ "MatRealPart_SeqDense" 1657 static PetscErrorCode MatRealPart_SeqDense(Mat A) 1658 { 1659 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1660 PetscInt i,nz = A->rmap->n*A->cmap->n; 1661 PetscScalar *aa = a->v; 1662 1663 PetscFunctionBegin; 1664 for (i=0; i<nz; i++) aa[i] = PetscRealPart(aa[i]); 1665 PetscFunctionReturn(0); 1666 } 1667 1668 #undef __FUNCT__ 1669 #define __FUNCT__ "MatImaginaryPart_SeqDense" 1670 static PetscErrorCode MatImaginaryPart_SeqDense(Mat A) 1671 { 1672 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1673 PetscInt i,nz = A->rmap->n*A->cmap->n; 1674 PetscScalar *aa = a->v; 1675 1676 PetscFunctionBegin; 1677 for (i=0; i<nz; i++) aa[i] = PetscImaginaryPart(aa[i]); 1678 PetscFunctionReturn(0); 1679 } 1680 1681 /* ----------------------------------------------------------------*/ 1682 #undef __FUNCT__ 1683 #define __FUNCT__ "MatMatMult_SeqDense_SeqDense" 1684 PetscErrorCode MatMatMult_SeqDense_SeqDense(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C) 1685 { 1686 PetscErrorCode ierr; 1687 1688 PetscFunctionBegin; 1689 if (scall == MAT_INITIAL_MATRIX) { 1690 ierr = MatMatMultSymbolic_SeqDense_SeqDense(A,B,fill,C);CHKERRQ(ierr); 1691 } 1692 ierr = MatMatMultNumeric_SeqDense_SeqDense(A,B,*C);CHKERRQ(ierr); 1693 PetscFunctionReturn(0); 1694 } 1695 1696 #undef __FUNCT__ 1697 #define __FUNCT__ "MatMatMultSymbolic_SeqDense_SeqDense" 1698 PetscErrorCode MatMatMultSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat *C) 1699 { 1700 PetscErrorCode ierr; 1701 PetscInt m=A->rmap->n,n=B->cmap->n; 1702 Mat Cmat; 1703 1704 PetscFunctionBegin; 1705 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); 1706 ierr = MatCreate(PETSC_COMM_SELF,&Cmat);CHKERRQ(ierr); 1707 ierr = MatSetSizes(Cmat,m,n,m,n);CHKERRQ(ierr); 1708 ierr = MatSetType(Cmat,MATSEQDENSE);CHKERRQ(ierr); 1709 ierr = MatSeqDenseSetPreallocation(Cmat,PETSC_NULL);CHKERRQ(ierr); 1710 1711 *C = Cmat; 1712 PetscFunctionReturn(0); 1713 } 1714 1715 #undef __FUNCT__ 1716 #define __FUNCT__ "MatMatMultNumeric_SeqDense_SeqDense" 1717 PetscErrorCode MatMatMultNumeric_SeqDense_SeqDense(Mat A,Mat B,Mat C) 1718 { 1719 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1720 Mat_SeqDense *b = (Mat_SeqDense*)B->data; 1721 Mat_SeqDense *c = (Mat_SeqDense*)C->data; 1722 PetscBLASInt m,n,k; 1723 PetscScalar _DOne=1.0,_DZero=0.0; 1724 PetscErrorCode ierr; 1725 1726 PetscFunctionBegin; 1727 ierr = PetscBLASIntCast(A->rmap->n,&m);CHKERRQ(ierr); 1728 ierr = PetscBLASIntCast(B->cmap->n,&n);CHKERRQ(ierr); 1729 ierr = PetscBLASIntCast(A->cmap->n,&k);CHKERRQ(ierr); 1730 BLASgemm_("N","N",&m,&n,&k,&_DOne,a->v,&a->lda,b->v,&b->lda,&_DZero,c->v,&c->lda); 1731 PetscFunctionReturn(0); 1732 } 1733 1734 #undef __FUNCT__ 1735 #define __FUNCT__ "MatTransposeMatMult_SeqDense_SeqDense" 1736 PetscErrorCode MatTransposeMatMult_SeqDense_SeqDense(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C) 1737 { 1738 PetscErrorCode ierr; 1739 1740 PetscFunctionBegin; 1741 if (scall == MAT_INITIAL_MATRIX) { 1742 ierr = MatTransposeMatMultSymbolic_SeqDense_SeqDense(A,B,fill,C);CHKERRQ(ierr); 1743 } 1744 ierr = MatTransposeMatMultNumeric_SeqDense_SeqDense(A,B,*C);CHKERRQ(ierr); 1745 PetscFunctionReturn(0); 1746 } 1747 1748 #undef __FUNCT__ 1749 #define __FUNCT__ "MatTransposeMatMultSymbolic_SeqDense_SeqDense" 1750 PetscErrorCode MatTransposeMatMultSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat *C) 1751 { 1752 PetscErrorCode ierr; 1753 PetscInt m=A->cmap->n,n=B->cmap->n; 1754 Mat Cmat; 1755 1756 PetscFunctionBegin; 1757 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); 1758 ierr = MatCreate(PETSC_COMM_SELF,&Cmat);CHKERRQ(ierr); 1759 ierr = MatSetSizes(Cmat,m,n,m,n);CHKERRQ(ierr); 1760 ierr = MatSetType(Cmat,MATSEQDENSE);CHKERRQ(ierr); 1761 ierr = MatSeqDenseSetPreallocation(Cmat,PETSC_NULL);CHKERRQ(ierr); 1762 1763 Cmat->assembled = PETSC_TRUE; 1764 1765 *C = Cmat; 1766 PetscFunctionReturn(0); 1767 } 1768 1769 #undef __FUNCT__ 1770 #define __FUNCT__ "MatTransposeMatMultNumeric_SeqDense_SeqDense" 1771 PetscErrorCode MatTransposeMatMultNumeric_SeqDense_SeqDense(Mat A,Mat B,Mat C) 1772 { 1773 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1774 Mat_SeqDense *b = (Mat_SeqDense*)B->data; 1775 Mat_SeqDense *c = (Mat_SeqDense*)C->data; 1776 PetscBLASInt m,n,k; 1777 PetscScalar _DOne=1.0,_DZero=0.0; 1778 PetscErrorCode ierr; 1779 1780 PetscFunctionBegin; 1781 ierr = PetscBLASIntCast(A->cmap->n,&m);CHKERRQ(ierr); 1782 ierr = PetscBLASIntCast(B->cmap->n,&n);CHKERRQ(ierr); 1783 ierr = PetscBLASIntCast(A->rmap->n,&k);CHKERRQ(ierr); 1784 /* 1785 Note the m and n arguments below are the number rows and columns of A', not A! 1786 */ 1787 BLASgemm_("T","N",&m,&n,&k,&_DOne,a->v,&a->lda,b->v,&b->lda,&_DZero,c->v,&c->lda); 1788 PetscFunctionReturn(0); 1789 } 1790 1791 #undef __FUNCT__ 1792 #define __FUNCT__ "MatGetRowMax_SeqDense" 1793 PetscErrorCode MatGetRowMax_SeqDense(Mat A,Vec v,PetscInt idx[]) 1794 { 1795 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1796 PetscErrorCode ierr; 1797 PetscInt i,j,m = A->rmap->n,n = A->cmap->n,p; 1798 PetscScalar *x; 1799 MatScalar *aa = a->v; 1800 1801 PetscFunctionBegin; 1802 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 1803 1804 ierr = VecSet(v,0.0);CHKERRQ(ierr); 1805 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 1806 ierr = VecGetLocalSize(v,&p);CHKERRQ(ierr); 1807 if (p != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 1808 for (i=0; i<m; i++) { 1809 x[i] = aa[i]; if (idx) idx[i] = 0; 1810 for (j=1; j<n; j++) { 1811 if (PetscRealPart(x[i]) < PetscRealPart(aa[i+m*j])) {x[i] = aa[i + m*j]; if (idx) idx[i] = j;} 1812 } 1813 } 1814 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 1815 PetscFunctionReturn(0); 1816 } 1817 1818 #undef __FUNCT__ 1819 #define __FUNCT__ "MatGetRowMaxAbs_SeqDense" 1820 PetscErrorCode MatGetRowMaxAbs_SeqDense(Mat A,Vec v,PetscInt idx[]) 1821 { 1822 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1823 PetscErrorCode ierr; 1824 PetscInt i,j,m = A->rmap->n,n = A->cmap->n,p; 1825 PetscScalar *x; 1826 PetscReal atmp; 1827 MatScalar *aa = a->v; 1828 1829 PetscFunctionBegin; 1830 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 1831 1832 ierr = VecSet(v,0.0);CHKERRQ(ierr); 1833 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 1834 ierr = VecGetLocalSize(v,&p);CHKERRQ(ierr); 1835 if (p != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 1836 for (i=0; i<m; i++) { 1837 x[i] = PetscAbsScalar(aa[i]); 1838 for (j=1; j<n; j++) { 1839 atmp = PetscAbsScalar(aa[i+m*j]); 1840 if (PetscAbsScalar(x[i]) < atmp) {x[i] = atmp; if (idx) idx[i] = j;} 1841 } 1842 } 1843 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 1844 PetscFunctionReturn(0); 1845 } 1846 1847 #undef __FUNCT__ 1848 #define __FUNCT__ "MatGetRowMin_SeqDense" 1849 PetscErrorCode MatGetRowMin_SeqDense(Mat A,Vec v,PetscInt idx[]) 1850 { 1851 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1852 PetscErrorCode ierr; 1853 PetscInt i,j,m = A->rmap->n,n = A->cmap->n,p; 1854 PetscScalar *x; 1855 MatScalar *aa = a->v; 1856 1857 PetscFunctionBegin; 1858 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 1859 1860 ierr = VecSet(v,0.0);CHKERRQ(ierr); 1861 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 1862 ierr = VecGetLocalSize(v,&p);CHKERRQ(ierr); 1863 if (p != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 1864 for (i=0; i<m; i++) { 1865 x[i] = aa[i]; if (idx) idx[i] = 0; 1866 for (j=1; j<n; j++) { 1867 if (PetscRealPart(x[i]) > PetscRealPart(aa[i+m*j])) {x[i] = aa[i + m*j]; if (idx) idx[i] = j;} 1868 } 1869 } 1870 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 1871 PetscFunctionReturn(0); 1872 } 1873 1874 #undef __FUNCT__ 1875 #define __FUNCT__ "MatGetColumnVector_SeqDense" 1876 PetscErrorCode MatGetColumnVector_SeqDense(Mat A,Vec v,PetscInt col) 1877 { 1878 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1879 PetscErrorCode ierr; 1880 PetscScalar *x; 1881 1882 PetscFunctionBegin; 1883 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 1884 1885 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 1886 ierr = PetscMemcpy(x,a->v+col*a->lda,A->rmap->n*sizeof(PetscScalar));CHKERRQ(ierr); 1887 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 1888 PetscFunctionReturn(0); 1889 } 1890 1891 1892 #undef __FUNCT__ 1893 #define __FUNCT__ "MatGetColumnNorms_SeqDense" 1894 PetscErrorCode MatGetColumnNorms_SeqDense(Mat A,NormType type,PetscReal *norms) 1895 { 1896 PetscErrorCode ierr; 1897 PetscInt i,j,m,n; 1898 PetscScalar *a; 1899 1900 PetscFunctionBegin; 1901 ierr = MatGetSize(A,&m,&n);CHKERRQ(ierr); 1902 ierr = PetscMemzero(norms,n*sizeof(PetscReal));CHKERRQ(ierr); 1903 ierr = MatDenseGetArray(A,&a);CHKERRQ(ierr); 1904 if (type == NORM_2) { 1905 for (i=0; i<n; i++) { 1906 for (j=0; j<m; j++) { 1907 norms[i] += PetscAbsScalar(a[j]*a[j]); 1908 } 1909 a += m; 1910 } 1911 } else if (type == NORM_1) { 1912 for (i=0; i<n; i++) { 1913 for (j=0; j<m; j++) { 1914 norms[i] += PetscAbsScalar(a[j]); 1915 } 1916 a += m; 1917 } 1918 } else if (type == NORM_INFINITY) { 1919 for (i=0; i<n; i++) { 1920 for (j=0; j<m; j++) { 1921 norms[i] = PetscMax(PetscAbsScalar(a[j]),norms[i]); 1922 } 1923 a += m; 1924 } 1925 } else SETERRQ(((PetscObject)A)->comm,PETSC_ERR_ARG_WRONG,"Unknown NormType"); 1926 ierr = MatDenseRestoreArray(A,&a);CHKERRQ(ierr); 1927 if (type == NORM_2) { 1928 for (i=0; i<n; i++) norms[i] = PetscSqrtReal(norms[i]); 1929 } 1930 PetscFunctionReturn(0); 1931 } 1932 1933 #undef __FUNCT__ 1934 #define __FUNCT__ "MatSetRandom_SeqDense" 1935 static PetscErrorCode MatSetRandom_SeqDense(Mat x,PetscRandom rctx) 1936 { 1937 PetscErrorCode ierr; 1938 PetscScalar *a; 1939 PetscInt m,n,i; 1940 1941 PetscFunctionBegin; 1942 ierr = MatGetSize(x,&m,&n);CHKERRQ(ierr); 1943 ierr = MatDenseGetArray(x,&a);CHKERRQ(ierr); 1944 for (i=0; i<m*n; i++) { 1945 ierr = PetscRandomGetValue(rctx,a+i);CHKERRQ(ierr); 1946 } 1947 ierr = MatDenseRestoreArray(x,&a);CHKERRQ(ierr); 1948 PetscFunctionReturn(0); 1949 } 1950 1951 1952 /* -------------------------------------------------------------------*/ 1953 static struct _MatOps MatOps_Values = { MatSetValues_SeqDense, 1954 MatGetRow_SeqDense, 1955 MatRestoreRow_SeqDense, 1956 MatMult_SeqDense, 1957 /* 4*/ MatMultAdd_SeqDense, 1958 MatMultTranspose_SeqDense, 1959 MatMultTransposeAdd_SeqDense, 1960 0, 1961 0, 1962 0, 1963 /* 10*/ 0, 1964 MatLUFactor_SeqDense, 1965 MatCholeskyFactor_SeqDense, 1966 MatSOR_SeqDense, 1967 MatTranspose_SeqDense, 1968 /* 15*/ MatGetInfo_SeqDense, 1969 MatEqual_SeqDense, 1970 MatGetDiagonal_SeqDense, 1971 MatDiagonalScale_SeqDense, 1972 MatNorm_SeqDense, 1973 /* 20*/ MatAssemblyBegin_SeqDense, 1974 MatAssemblyEnd_SeqDense, 1975 MatSetOption_SeqDense, 1976 MatZeroEntries_SeqDense, 1977 /* 24*/ MatZeroRows_SeqDense, 1978 0, 1979 0, 1980 0, 1981 0, 1982 /* 29*/ MatSetUp_SeqDense, 1983 0, 1984 0, 1985 0, 1986 0, 1987 /* 34*/ MatDuplicate_SeqDense, 1988 0, 1989 0, 1990 0, 1991 0, 1992 /* 39*/ MatAXPY_SeqDense, 1993 MatGetSubMatrices_SeqDense, 1994 0, 1995 MatGetValues_SeqDense, 1996 MatCopy_SeqDense, 1997 /* 44*/ MatGetRowMax_SeqDense, 1998 MatScale_SeqDense, 1999 0, 2000 0, 2001 0, 2002 /* 49*/ MatSetRandom_SeqDense, 2003 0, 2004 0, 2005 0, 2006 0, 2007 /* 54*/ 0, 2008 0, 2009 0, 2010 0, 2011 0, 2012 /* 59*/ 0, 2013 MatDestroy_SeqDense, 2014 MatView_SeqDense, 2015 0, 2016 0, 2017 /* 64*/ 0, 2018 0, 2019 0, 2020 0, 2021 0, 2022 /* 69*/ MatGetRowMaxAbs_SeqDense, 2023 0, 2024 0, 2025 0, 2026 0, 2027 /* 74*/ 0, 2028 0, 2029 0, 2030 0, 2031 0, 2032 /* 79*/ 0, 2033 0, 2034 0, 2035 0, 2036 /* 83*/ MatLoad_SeqDense, 2037 0, 2038 MatIsHermitian_SeqDense, 2039 0, 2040 0, 2041 0, 2042 /* 89*/ MatMatMult_SeqDense_SeqDense, 2043 MatMatMultSymbolic_SeqDense_SeqDense, 2044 MatMatMultNumeric_SeqDense_SeqDense, 2045 0, 2046 0, 2047 /* 94*/ 0, 2048 0, 2049 0, 2050 0, 2051 0, 2052 /* 99*/ 0, 2053 0, 2054 0, 2055 MatConjugate_SeqDense, 2056 0, 2057 /*104*/ 0, 2058 MatRealPart_SeqDense, 2059 MatImaginaryPart_SeqDense, 2060 0, 2061 0, 2062 /*109*/ MatMatSolve_SeqDense, 2063 0, 2064 MatGetRowMin_SeqDense, 2065 MatGetColumnVector_SeqDense, 2066 0, 2067 /*114*/ 0, 2068 0, 2069 0, 2070 0, 2071 0, 2072 /*119*/ 0, 2073 0, 2074 0, 2075 0, 2076 0, 2077 /*124*/ 0, 2078 MatGetColumnNorms_SeqDense, 2079 0, 2080 0, 2081 0, 2082 /*129*/ 0, 2083 MatTransposeMatMult_SeqDense_SeqDense, 2084 MatTransposeMatMultSymbolic_SeqDense_SeqDense, 2085 MatTransposeMatMultNumeric_SeqDense_SeqDense, 2086 }; 2087 2088 #undef __FUNCT__ 2089 #define __FUNCT__ "MatCreateSeqDense" 2090 /*@C 2091 MatCreateSeqDense - Creates a sequential dense matrix that 2092 is stored in column major order (the usual Fortran 77 manner). Many 2093 of the matrix operations use the BLAS and LAPACK routines. 2094 2095 Collective on MPI_Comm 2096 2097 Input Parameters: 2098 + comm - MPI communicator, set to PETSC_COMM_SELF 2099 . m - number of rows 2100 . n - number of columns 2101 - data - optional location of matrix data in column major order. Set data=PETSC_NULL for PETSc 2102 to control all matrix memory allocation. 2103 2104 Output Parameter: 2105 . A - the matrix 2106 2107 Notes: 2108 The data input variable is intended primarily for Fortran programmers 2109 who wish to allocate their own matrix memory space. Most users should 2110 set data=PETSC_NULL. 2111 2112 Level: intermediate 2113 2114 .keywords: dense, matrix, LAPACK, BLAS 2115 2116 .seealso: MatCreate(), MatCreateDense(), MatSetValues() 2117 @*/ 2118 PetscErrorCode MatCreateSeqDense(MPI_Comm comm,PetscInt m,PetscInt n,PetscScalar *data,Mat *A) 2119 { 2120 PetscErrorCode ierr; 2121 2122 PetscFunctionBegin; 2123 ierr = MatCreate(comm,A);CHKERRQ(ierr); 2124 ierr = MatSetSizes(*A,m,n,m,n);CHKERRQ(ierr); 2125 ierr = MatSetType(*A,MATSEQDENSE);CHKERRQ(ierr); 2126 ierr = MatSeqDenseSetPreallocation(*A,data);CHKERRQ(ierr); 2127 PetscFunctionReturn(0); 2128 } 2129 2130 #undef __FUNCT__ 2131 #define __FUNCT__ "MatSeqDenseSetPreallocation" 2132 /*@C 2133 MatSeqDenseSetPreallocation - Sets the array used for storing the matrix elements 2134 2135 Collective on MPI_Comm 2136 2137 Input Parameters: 2138 + A - the matrix 2139 - data - the array (or PETSC_NULL) 2140 2141 Notes: 2142 The data input variable is intended primarily for Fortran programmers 2143 who wish to allocate their own matrix memory space. Most users should 2144 need not call this routine. 2145 2146 Level: intermediate 2147 2148 .keywords: dense, matrix, LAPACK, BLAS 2149 2150 .seealso: MatCreate(), MatCreateDense(), MatSetValues(), MatSeqDenseSetLDA() 2151 2152 @*/ 2153 PetscErrorCode MatSeqDenseSetPreallocation(Mat B,PetscScalar data[]) 2154 { 2155 PetscErrorCode ierr; 2156 2157 PetscFunctionBegin; 2158 ierr = PetscTryMethod(B,"MatSeqDenseSetPreallocation_C",(Mat,PetscScalar[]),(B,data));CHKERRQ(ierr); 2159 PetscFunctionReturn(0); 2160 } 2161 2162 EXTERN_C_BEGIN 2163 #undef __FUNCT__ 2164 #define __FUNCT__ "MatSeqDenseSetPreallocation_SeqDense" 2165 PetscErrorCode MatSeqDenseSetPreallocation_SeqDense(Mat B,PetscScalar *data) 2166 { 2167 Mat_SeqDense *b; 2168 PetscErrorCode ierr; 2169 2170 PetscFunctionBegin; 2171 B->preallocated = PETSC_TRUE; 2172 2173 ierr = PetscLayoutSetUp(B->rmap);CHKERRQ(ierr); 2174 ierr = PetscLayoutSetUp(B->cmap);CHKERRQ(ierr); 2175 2176 b = (Mat_SeqDense*)B->data; 2177 b->Mmax = B->rmap->n; 2178 b->Nmax = B->cmap->n; 2179 if (b->lda <= 0 || b->changelda) b->lda = B->rmap->n; 2180 2181 if (!data) { /* petsc-allocated storage */ 2182 if (!b->user_alloc) { ierr = PetscFree(b->v);CHKERRQ(ierr); } 2183 ierr = PetscMalloc(b->lda*b->Nmax*sizeof(PetscScalar),&b->v);CHKERRQ(ierr); 2184 ierr = PetscMemzero(b->v,b->lda*b->Nmax*sizeof(PetscScalar));CHKERRQ(ierr); 2185 ierr = PetscLogObjectMemory(B,b->lda*b->Nmax*sizeof(PetscScalar));CHKERRQ(ierr); 2186 2187 b->user_alloc = PETSC_FALSE; 2188 } else { /* user-allocated storage */ 2189 if (!b->user_alloc) { ierr = PetscFree(b->v);CHKERRQ(ierr); } 2190 b->v = data; 2191 b->user_alloc = PETSC_TRUE; 2192 } 2193 B->assembled = PETSC_TRUE; 2194 PetscFunctionReturn(0); 2195 } 2196 EXTERN_C_END 2197 2198 #undef __FUNCT__ 2199 #define __FUNCT__ "MatSeqDenseSetLDA" 2200 /*@C 2201 MatSeqDenseSetLDA - Declare the leading dimension of the user-provided array 2202 2203 Input parameter: 2204 + A - the matrix 2205 - lda - the leading dimension 2206 2207 Notes: 2208 This routine is to be used in conjunction with MatSeqDenseSetPreallocation(); 2209 it asserts that the preallocation has a leading dimension (the LDA parameter 2210 of Blas and Lapack fame) larger than M, the first dimension of the matrix. 2211 2212 Level: intermediate 2213 2214 .keywords: dense, matrix, LAPACK, BLAS 2215 2216 .seealso: MatCreate(), MatCreateSeqDense(), MatSeqDenseSetPreallocation(), MatSetMaximumSize() 2217 2218 @*/ 2219 PetscErrorCode MatSeqDenseSetLDA(Mat B,PetscInt lda) 2220 { 2221 Mat_SeqDense *b = (Mat_SeqDense*)B->data; 2222 2223 PetscFunctionBegin; 2224 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); 2225 b->lda = lda; 2226 b->changelda = PETSC_FALSE; 2227 b->Mmax = PetscMax(b->Mmax,lda); 2228 PetscFunctionReturn(0); 2229 } 2230 2231 /*MC 2232 MATSEQDENSE - MATSEQDENSE = "seqdense" - A matrix type to be used for sequential dense matrices. 2233 2234 Options Database Keys: 2235 . -mat_type seqdense - sets the matrix type to "seqdense" during a call to MatSetFromOptions() 2236 2237 Level: beginner 2238 2239 .seealso: MatCreateSeqDense() 2240 2241 M*/ 2242 2243 EXTERN_C_BEGIN 2244 #undef __FUNCT__ 2245 #define __FUNCT__ "MatCreate_SeqDense" 2246 PetscErrorCode MatCreate_SeqDense(Mat B) 2247 { 2248 Mat_SeqDense *b; 2249 PetscErrorCode ierr; 2250 PetscMPIInt size; 2251 2252 PetscFunctionBegin; 2253 ierr = MPI_Comm_size(((PetscObject)B)->comm,&size);CHKERRQ(ierr); 2254 if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Comm must be of size 1"); 2255 2256 ierr = PetscNewLog(B,Mat_SeqDense,&b);CHKERRQ(ierr); 2257 ierr = PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 2258 B->data = (void*)b; 2259 2260 b->pivots = 0; 2261 b->roworiented = PETSC_TRUE; 2262 b->v = 0; 2263 b->changelda = PETSC_FALSE; 2264 2265 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatDenseGetArray_C","MatDenseGetArray_SeqDense",MatDenseGetArray_SeqDense);CHKERRQ(ierr); 2266 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatDenseRestoreArray_C","MatDenseRestoreArray_SeqDense",MatDenseRestoreArray_SeqDense);CHKERRQ(ierr); 2267 2268 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqdense_seqaij_C","MatConvert_SeqDense_SeqAIJ",MatConvert_SeqDense_SeqAIJ);CHKERRQ(ierr); 2269 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_petsc_C", 2270 "MatGetFactor_seqdense_petsc", 2271 MatGetFactor_seqdense_petsc);CHKERRQ(ierr); 2272 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqDenseSetPreallocation_C", 2273 "MatSeqDenseSetPreallocation_SeqDense", 2274 MatSeqDenseSetPreallocation_SeqDense);CHKERRQ(ierr); 2275 2276 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatMatMult_seqaij_seqdense_C", 2277 "MatMatMult_SeqAIJ_SeqDense", 2278 MatMatMult_SeqAIJ_SeqDense);CHKERRQ(ierr); 2279 2280 2281 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatMatMultSymbolic_seqaij_seqdense_C", 2282 "MatMatMultSymbolic_SeqAIJ_SeqDense", 2283 MatMatMultSymbolic_SeqAIJ_SeqDense);CHKERRQ(ierr); 2284 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatMatMultNumeric_seqaij_seqdense_C", 2285 "MatMatMultNumeric_SeqAIJ_SeqDense", 2286 MatMatMultNumeric_SeqAIJ_SeqDense);CHKERRQ(ierr); 2287 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQDENSE);CHKERRQ(ierr); 2288 PetscFunctionReturn(0); 2289 } 2290 EXTERN_C_END 2291