Actual source code: ex243.c

  1: static char help[] = "Test conversion of ScaLAPACK matrices.\n\n";

  3: #include <petscmat.h>

  5: int main(int argc, char** argv)
  6: {
  7:   Mat            A,A_scalapack;
  8:   PetscInt       i,j,M=10,N=5,nloc,mloc,nrows,ncols;
 10:   PetscMPIInt    rank,size;
 11:   IS             isrows,iscols;
 12:   const PetscInt *rows,*cols;
 13:   PetscScalar    *v;
 14:   MatType        type;
 15:   PetscBool      isDense,isAIJ,flg;

 17:   PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
 18:   MPI_Comm_size(PETSC_COMM_WORLD,&size);
 19:   MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
 20:   PetscOptionsGetInt(NULL,NULL,"-M",&M,NULL);
 21:   PetscOptionsGetInt(NULL,NULL,"-N",&N,NULL);

 23:   /* Create a matrix */
 24:   MatCreate(PETSC_COMM_WORLD, &A);
 25:   mloc = PETSC_DECIDE;
 26:   PetscSplitOwnershipEqual(PETSC_COMM_WORLD,&mloc,&M);
 27:   nloc = PETSC_DECIDE;
 28:   PetscSplitOwnershipEqual(PETSC_COMM_WORLD,&nloc,&N);
 29:   MatSetSizes(A,mloc,nloc,M,N);
 30:   MatSetType(A,MATDENSE);
 31:   MatSetFromOptions(A);
 32:   MatSetUp(A);

 34:   /* Set local matrix entries */
 35:   MatGetOwnershipIS(A,&isrows,&iscols);
 36:   ISGetLocalSize(isrows,&nrows);
 37:   ISGetIndices(isrows,&rows);
 38:   ISGetLocalSize(iscols,&ncols);
 39:   ISGetIndices(iscols,&cols);
 40:   PetscMalloc1(nrows*ncols,&v);

 42:   for (i=0; i<nrows; i++) {
 43:     for (j=0; j<ncols; j++) {
 44:       if (size == 1) {
 45:         v[i*ncols+j] = (PetscScalar)(i+j);
 46:       } else {
 47:         v[i*ncols+j] = (PetscScalar)rank+j*0.1;
 48:       }
 49:     }
 50:   }
 51:   MatSetValues(A,nrows,rows,ncols,cols,v,INSERT_VALUES);
 52:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
 53:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);

 55:   /* Test MatSetValues() by converting A to A_scalapack */
 56:   MatGetType(A,&type);
 57:   if (size == 1) {
 58:     PetscObjectTypeCompare((PetscObject)A,MATSEQDENSE,&isDense);
 59:     PetscObjectTypeCompare((PetscObject)A,MATSEQAIJ,&isAIJ);
 60:   } else {
 61:     PetscObjectTypeCompare((PetscObject)A,MATMPIDENSE,&isDense);
 62:     PetscObjectTypeCompare((PetscObject)A,MATMPIAIJ,&isAIJ);
 63:   }

 65:   if (isDense || isAIJ) {
 66:     Mat Aexplicit;
 67:     MatConvert(A,MATSCALAPACK,MAT_INITIAL_MATRIX,&A_scalapack);
 68:     MatComputeOperator(A_scalapack,isAIJ?MATAIJ:MATDENSE,&Aexplicit);
 69:     MatMultEqual(Aexplicit,A_scalapack,5,&flg);
 70:     if (!flg) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_PLIB,"Aexplicit != A_scalapack.");
 71:     MatDestroy(&Aexplicit);

 73:     /* Test MAT_REUSE_MATRIX which is only supported for inplace conversion */
 74:     MatConvert(A,MATSCALAPACK,MAT_INPLACE_MATRIX,&A);
 75:     MatMultEqual(A_scalapack,A,5,&flg);
 76:     if (!flg) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_PLIB,"A_scalapack != A.");
 77:     MatDestroy(&A_scalapack);
 78:   }

 80:   ISRestoreIndices(isrows,&rows);
 81:   ISRestoreIndices(iscols,&cols);
 82:   ISDestroy(&isrows);
 83:   ISDestroy(&iscols);
 84:   PetscFree(v);
 85:   MatDestroy(&A);
 86:   PetscFinalize();
 87:   return ierr;
 88: }

 90: /*TEST

 92:    build:
 93:       requires: scalapack

 95:    test:
 96:       nsize: 6

 98:    test:
 99:       suffix: 2
100:       nsize: 6
101:       args: -mat_type aij
102:       output_file: output/ex243_1.out

104:    test:
105:       suffix: 3
106:       nsize: 6
107:       args: -mat_type scalapack
108:       output_file: output/ex243_1.out

110: TEST*/