[petsc-users] 3-dimension to matrix

[Joon hee Choi]

Hi Jed Brown,

Thank you for your fast reply. Your last comments looks like my first code. The full size of the matrix is (X, Y*Z). Also, I used SEQAIJ as matrix type and (X, Y) as block size. Also, I implemented SeqAIJpreallocating with nnz. Nevertheless, it was very slow. The Matrix-Set-Up part of the first code is as follows:

        sort(tups.begin(), tups.end());
        MatCreate(PETSC_COMM_SELF, &A);
        MatSetType(A, MATSEQAIJ);
        MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, X, Y*Z);
        MatSetBlockSizes(A, X, Y);
        MatSeqAIJSetPreallocation(A, PETSC_DEFAULT, nnz);

        sz = tups.size();
        for (i=0; i<sz; i++) {
                x = std::tr1::get<0>(tups[i]);
                y = std::tr1::get<2>(tups[i]) + std::tr1::get<1>(tups[i])*Y;
                val = std::tr1::get<3>(tups[i]);
                MatSetValues(A, 1, &x, 1, &y, &val, INSERT_VALUES);
        }
        MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
        MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);


I used tuple (x, z, y, value) and vector of c++. I didn't get any errors from this code. However, it took about 9 minutes in this part. Please let me know what to change my code. Thank you.

Joon

----- Original Message -----
From: "Jed Brown" <jedbrown at mcs.anl.gov>
To: "Joon hee Choi" <choi240 at purdue.edu>, petsc-users at mcs.anl.gov
Sent: Friday, May 17, 2013 8:33:18 AM
Subject: Re: [petsc-users] 3-dimension to matrix

Joon hee Choi <choi240 at purdue.edu> writes:

> Hi all,
>
> I am trying to set up a matrix with 3d coordinates(x,y,z) and
> value. If X=max(x), Y=max(y) and Z=max(z), then the matrix is as
> follows: A = [ A1 A2 ... AZ ] (each A* is X-by-Y block matrix and A
> consists of Z block matrices)

Are these decoupled diagonal blocks or a sequence of blocks in the row.
That is, is the full size of the matrix (X, Y*Z) or (X*Z, Y*Z).  From
the code below, it looks like you intend for it to be the former.  What
parallel distribution do you intend for the input vector (of dimension
Y*Z) and the output vector (of dimension X)?

> When I tried to set up the matrix using MatSetValues with each
> coordinates and value, I could get the matrix, but it took 10
> minutes. 

This is almost certainly due to not preallocating.

> So I tried to set up the matrix again using
> MatSetLocalToGlobalMapping(Mat x, ISLocalToGlobalMapping rmapping,
> ISLocalToGlobalMapping cmapping), but ISLocalToGloalMapping was minus
> because my dataset is too big (X:26M, Y:26M, Z:48M, nonzero:144M) and
> I got errors. Is there the way to set up the matrix quickly in petsc?
> If so, please let me know the way. And my code using
> MatSetLocalToGlobalMapping following:
>
>
> ISLocalToGlobalMappingCreate(PETSC_COMM_SELF, 1, &zr, PETSC_OWN_POINTER, &irow);

This has only one entry.

> ISLocalToGlobalMappingUnBlock(irow, X, &orow);

How did you set the block size of the matrix?

>
> for (z=0; z<Z; i++)
> {
>     ISLocalToGlobalMappingCreate(PETSC_COMM_SELF, 1, &z, PETSC_OWN_POINTER, &icol);

This mapping also has only one entry.  It is not intended for you to
change the local-to-global mapping on a matrix.

>     MatSetLocalToGlobalMappingBlock(A, irow, icol);
>     ISLocalToGlobalMappingUnBlock(icol, Y, &ocol);
>     MatSetLocalToGlobalMapping(A, orow, ocol);
>
>     while(!mapX[z].empty())
>     {
>        x = mapX[z].front();
>        y = mapY[z].front();

So instead of changing it, just add the offset to y here.

>        val = mapVal[z].front();
>        MatSetValuesLocal(A, 1, &x, 1, &y, &val, INSERT_VALUES);
>        mapX[z].pop_front();
>        mapY[z].pop_front();
>        mapVal[z].pop_front();
>     }
> }
>
> Thank you,
>
> Joon