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*/