No I also printed those out and checked them. They return 0->nGlobalAtoms. The globalID-1 is because the other library which is assigning these ID's (LAMMPS) doesn’t use 0-based indexing for the global indices for whatever reason.
Best regards, Jacob Faibussowitsch (Jacob Fai - booss - oh - vitch) Cell: (312) 694-3391 > On Apr 8, 2020, at 5:14 PM, Jed Brown <[email protected]> wrote: > > Sure, but could you have computed them as a negative number here? > > Jacob Faibussowitsch <[email protected]> writes: > >> Supposed to be the block row and column indices. So they are the atom global >> IDs. >> >> Best regards, >> >> Jacob Faibussowitsch >> (Jacob Fai - booss - oh - vitch) >> Cell: (312) 694-3391 >> >>> On Apr 8, 2020, at 4:29 PM, Jed Brown <[email protected]> wrote: >>> >>> What are idxm[0] and idxn[0] in this code? >>> >>> Jacob Faibussowitsch <[email protected]> writes: >>> >>>> 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
