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