Hello All,
I am using MATSBAIJ to make a symmetric hamiltonian matrix, but for some reason
MatSetValuesBlocked isn’t inserting anything.
Setup:
I want a 4n x 4n symmetric matrix, where each MatSetValuesBlocked inserts a 4x4
sub-matrix of values into the global matrix, positioned by the global ID of two
atoms that I am extracting. So the 4x4 interaction sub matrix between atom 1
and atom 4 will be in rows 0-4 and cols 16-20.
My non-working implementation:
PetscInt bs = 4;
ierr = MatCreateSBAIJ(comm, bs, bs*nLocalAtoms, bs*nLocalAtoms,
PETSC_DETERMINE, PETSC_DETERMINE, 1, NULL, 1, NULL, H);CHKERRQ(ierr);
ierr = MatSetOption(*H, MAT_SYMMETRIC, PETSC_TRUE);CHKERRQ(ierr);
ierr = MatSetUp(*H);CHKERRQ(ierr);
.
.
.
START LOOP {
Calculations...
const PetscScalar valmat[16] = {Es[0], Es[1], Es[2], Es[3], Ex[0], Ex[1],
Ex[2], Ex[3], Ey[0], Ey[1], Ey[2], Ey[3], Ez[0], Ez[1], Ez[2], Ez[3]};
PetscInt idxm[1], idxn[1];
PetscInt m = 1, n = 1;
if (globalIDs[i] < globalIDs[neighidx[0]]-1) {
// Insert in upper triangular
idxm[0] = globalIDs[i]-1;
idxn[0] = globalIDs[neighidx[0]]-1;
} else {
// swap so we insert on upper triangular
idxn[0] = globalIDs[i]-1;
idxm[0] = globalIDs[neighidx[0]]-1;
}
ierr = MatSetValuesBlocked(*H, m, (const PetscInt *) idxm, n, (const PetscInt
*) idxn, valmat, INSERT_VALUES);CHKERRQ(ierr);
} END LOOP
ierr = MatAssemblyBegin(*H, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
ierr = MatAssemblyEnd(*H, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
RETURN H
Output from -mat_view:
Mat Object: 2 MPI processes
type: mpisbaij
row 0:
row 1:
row 2:
row 3:
row 4:
row 5:
row 6:
row 7:
row 8:
row 9:
row 10:
row 11:
row 12:
row 13:
row 14:
row 15:
I also checked -info :mat: and it also reports X unneeded storage space but
more importantly, 0 used. I am not sure what is going wrong here.
Best regards,
Jacob Faibussowitsch
(Jacob Fai - booss - oh - vitch)
Cell: (312) 694-3391
