On Thu, Oct 10, 2013 at 11:50 AM, Christophe Ortiz < [email protected]> wrote:
> > On Wed, Oct 9, 2013 at 10:51 PM, Barry Smith <[email protected]> wrote: > >> >> On Oct 9, 2013, at 2:04 PM, Christophe Ortiz <[email protected]> >> wrote: >> >> > >> > On Wed, Oct 9, 2013 at 3:13 PM, Christophe Ortiz < >> [email protected]> wrote: >> > >> > > For instance, if we take the example from petsc (MatSetValuesBlocked >> page), >> > > what should I do to pass the block of values 3, 4, 7 and 8 with >> > > MatSetValuesBlocked ? >> > > >> > > 1 2 | 3 4 >> > > 5 6 | 7 8 >> > > - - - | - - - >> > > 9 10 | 11 12 >> > > 13 14 | 15 16 >> > >> > Block row 0, block column 1. (Count from 0.) >> > >> > >> > >> > And if with the previous case I would like to pass the values 1, 2, 3, >> 4, 5, 6, 7 and 8 (the two first rows) what is the matrix val[ ] that I >> should pass ? Should it be a matrix val[2][4] ? >> >> It takes a one dimensional array as input; by default the values >> should be ordered by row so you would pass in the numerical values 1, 2, 3, >> 4,5, 6,7,8 in that order. If you wish to have the values ordered by column >> then you would pass 1 5 2 6 3 7 4 8 after calling >> MatSetOptions(mat,MAT_ROW_ORIENTED,PETSC_FALSE). >> >> Note that the ordering you use on the numerical values is the same for >> both MatSetValues() and MatSetValuesBlocked(). >> > > > I am facing some problems with this when filling in the Jacobian for a > problem of 2 diffusion equations (dof=2). See ex25.c. > > If I'm right, then for node i, there should be 3 blocks of 2x2: > > X 0 | Y 0 | X 0 > 0 X | 0 Y | 0 X > > Is it correct ? > > First I tried what you said, to pass a one dimensional array: > row = i; > cols[0] = i-1; cols[1] = i; cols[2] = i+1; > > j = 0; > vallin[j] = X; vallin[j++] = 0; vallin[j++] = Y; vallin[j++] = 0; > vallin[j++] = X; vallin[j++] = 0; > > vallin[j++] = 0; vallin[j++] = X; vallin[j++] = 0; vallin[j++] = Y; > vallin[j++] = 0; vallin[j++] = X; > > ierr = MatSetValuesBlocked(*Jpre,1,&row,3,cols,&vallin[0],INSERT_VALUES); > > I hit the following message: > > [0]PETSC ERROR: --------------------- Error Message > ------------------------------------ > [0]PETSC ERROR: Floating point exception! > [0]PETSC ERROR: Vec entry at local location 0 is not-a-number or infinite > at end of function: Parameter number 3! > > > Then, I tried a 2-dimensional array. I saw in some examples of > ts/tutorials that it is also possible to pass a multi-dimensional arrays. > After all, in memory the rows are stored contiguously, so I guess it is > equivalent to a one-dimensional array. > > j = 0; k = 0; > valdiff[j][k] = X; valdiff[j][k++] = 0; > valdiff[j][k++] = Y; valdiff[j][k++] = 0; > valdiff[j][k++] = X; valdiff[j][k++] = 0; > > j= 1; k = 0; > valdiff[j][k] = 0; valdiff[j][k++] = X; > valdiff[j][k++] = 0; valdiff[j][k++] = Y; > valdiff[j][k++] = 0; valdiff[j][k++] = X; > > ierr = > MatSetValuesBlocked(*Jpre,1,&row,3,cols,&valdiff[0][0],INSERT_VALUES); > > I get the same error message. > What did I do wrong ? > > In ex25.c they use a 3-dimensional array: > > const PetscScalar vals[2][3][2] = {{{dxxL,0},{a *hx+dxx0,0},{dxxR,0}}, > > {{0,dxxL},{0,a*hx+dxx0},{0,dxxR}}}; > > and &vals[0][0][0] is passed in MatSetValuesBlocked (....,...., > &vals[0][0][0],...), and it works, though I don't understand how. > > > > I found my mistake ! In the indices of the arrays I should have written ++j and not j++. A good printf of each value helped me to find out...:-). Now it works. Christophe
