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