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