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
