What are idxm[0] and idxn[0] in this code? Jacob Faibussowitsch <jacob....@gmail.com> 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