#include<stdio.h>
#include<math.h>
#include<stdlib.h>
#include "petscvec.h" /* This enables us to use vectors. */
#include "petscmat.h" /* This enables us to use Matrices. It includes the petscvec header file*/
#include "petscksp.h" /* Now we can solve linear systems. Solvers used are KSP. */
int main(int argc, char **argv)
{
/* Declaration of Matrix A and some vectors x*/
Vec b;
int tempi,tempj;
Mat A;
FILE *fp;
MPI_Comm comm;
PetscInt n=2683,m=1274,index,lenB=39051;
PetscScalar scalar,rhs[n],B[lenB][3];
PetscErrorCode ierr;
PetscInt i,j;
comm = MPI_COMM_SELF;
/* This part is needed for the help flag supplied at run-time*/
ierr = PetscInitialize(&argc,&argv,(char*)0,help);CHKERRQ(ierr);
ierr = PetscOptionsGetInt(PETSC_NULL,"-n",&n,PETSC_NULL);CHKERRQ(ierr);
/* Use options from the terminal to create a vector that is type shared or mpi. */
ierr = VecCreate(PETSC_COMM_WORLD,&b);CHKERRQ(ierr); /* Vector creation */
ierr = VecSetSizes(b,PETSC_DECIDE,n);CHKERRQ(ierr); /* Setting the vector size */
ierr = VecSetFromOptions(b);CHKERRQ(ierr); /* Setting the vector type (shared, mpi etc) */
/* The second argument is a PETSc scalar value.*/
ierr = VecSet(b,0);CHKERRQ(ierr);
/* Reading in the RHS vector. */
/* Reading in the matrix from the file matrix.txt */
fp=fopen("b1.mat","r");
if (fp==NULL)
{
fprintf(stderr, "Cannot open file");
exit(1);
}
for (i = 0; i < n; i++)
{
if (fscanf(fp,"%lf", &rhs[i]) != 1)
{
fprintf(stderr, "Failed to read rhs vector[%d]\n", i);
exit(1);
}
}
index=0;
/*Putting x into final form */
for (i=0; i<n; ++i)
{
ierr= VecSetValues(b,1,&index,&rhs[i],INSERT_VALUES);CHKERRQ(ierr); /* One insertion per step. Thats what the 1 in second argument stands for */
index=index+1;
} /* The third and fourth arguments are addresses. The fifth argument is IORA */
/* Assembling the vector. */
ierr= VecAssemblyBegin(b);CHKERRQ(ierr);
ierr=VecAssemblyEnd(b);CHKERRQ(ierr);
/* Viewing the changed vector. */
ierr=PetscPrintf(PETSC_COMM_WORLD,"Vector b:\n");CHKERRQ(ierr);
ierr=VecView(b,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
/* ------------------------------------------------------------------------------------------------------------------------------------------------------ */
ierr=MatCreate(PETSC_COMM_WORLD,&A);CHKERRQ(ierr);
ierr = MatSetSizes(A,n,m,n,m);CHKERRQ(ierr); /* Setting the matrix size */
ierr = MatSetFromOptions(A);CHKERRQ(ierr); /* Setting the matrix type (shared, mpi etc) */
/* Reading in the matrix from the file matrix.txt */
fp=fopen("N1.mat","r");
if (fp==NULL)
{
fprintf(stderr, "Cannot open file");
exit(1);
}
for (i = 0; i < lenB; i++)
{
for (j = 0; j < 3; j++)
{
if (fscanf(fp,"%lf", &B[i][j]) != 1)
{
fprintf(stderr, "Failed to read matrix[%d][%d]\n", i, j);
exit(1);
}
}
}
/* Initializing the matrix A to zero matrix */
for(i=0;i<n;++i)
{
for(j=0;j<m;++j)
{
scalar=0.0;
ierr=MatSetValues(A,1,&i,1,&j,&scalar,INSERT_VALUES);CHKERRQ(ierr);
}
}
ierr=MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
ierr=MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
/* Inserting the non-zero elements. */
for(i=0;i<lenB;++i)
{
tempi=(int)B[i][0]-1;
tempj=(int)B[i][1]-1;
ierr=MatSetValues(A,1,&tempi,1,&tempj,&B[i][2],INSERT_VALUES);CHKERRQ(ierr);
}
ierr=MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
ierr=MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
//ierr=MatView(A,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
/*
Destroy any objects created.
*/
ierr=VecDestroy(b);CHKERRQ(ierr);
ierr=MatDestroy(A);CHKERRQ(ierr);
ierr=PetscFinalize();CHKERRQ(ierr);
return 0;
}