[petsc-users] Using KSPLSQR without explicitly computing U(transpose)U
Gaurish Telang
gaurish108 at gmail.com
Sun Jan 23 00:39:18 CST 2011
Hi,
I have been trying to solve an overdetermined system of linear equations on
PETSc using KSPLSQR.
Now for my minimization problem |Ux-b| where U has dimension 2683x1274, U
has full rank but is badly conditioned. This makes the algorithm LSQR an
ideal algorithm for doing Least sqaures here, (at least according to the
paper by Paige and Saunders.)
I have been able to solve the normal equations ,
U(transpose)Ux=U(transpose)b successfully with KSPLSQR, in the usual way one
solves linear systems with SQUARE matrices.
But for this I had to compute U(transpose)U explicitly. The algorithm
mentioned in the paper does NOT involve an explicit computation of
U(transpose)U.
Is there a way to avoid the explicit and expensive computation of
U(transpose)U and use KSPLSQR? The details on this page
http://www.mcs.anl.gov/petsc/petsc-2/snapshots/petsc-current/docs/manualpages/KSP/KSPLSQR.htmldid
not really help.
I have also experimented with the suggesstions in the thread
http://lists.mcs.anl.gov/pipermail/petsc-users/2010-August/006784.html but I
kept getting errors saying there is a dimension mismatch
I use KSPSetOperators as ierr =
KSPSetOperators(ksp,Prod,Prod,SAME_PRECONDITIONER);CHKERRQ(ierr); where
Prod=U(Transpose)U calculated in the previous statements.
Anyway, I am providing my current working code underneath just for
reference. U and b are provided through files.
I am guessing the answer lies in the way one uses KSPSetOperators but I am
not sure.
Regards,
Gaurish Telang
%--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
static char help[] = " Doing least squares \n\
-f <input_file> : file to load \n\n";
/*
Include "petscmat.h" so that we can use matrices.
automatically includes:
petscsys.h - base PETSc routines petscvec.h - vectors
petscmat.h - matrices
petscis.h - index sets petscviewer.h -
viewers
*/
#include "petscmat.h"
#include "petscvec.h"
#include "petscksp.h" /* For the iterative solvers */
#include <stdlib.h>
#include <stdio.h>
#undef __FUNCT__
#define __FUNCT__ "main"
int main(int argc,char **args)
{
Mat U,U_Transpose,Prod; /* matrix */
PetscViewer fd; /* viewer */
char file[PETSC_MAX_PATH_LEN]; /* input file name */
PetscErrorCode ierr;
PetscTruth flg;
Vec b,x,new_rhs;
PetscInt i,n,m,bsize=2683,index; /* bsize indicates the text
file size of bx.mat */
PetscScalar *xx;
PetscScalar rhs[bsize];
FILE *fp;
KSP ksp;
PC pc;
PetscMPIInt size;
PetscInt num_iters;
PetscReal rnorm;
KSPConvergedReason reason;
PetscInitialize(&argc,&args,(char *)0,help);
ierr =
PetscOptionsGetString(PETSC_NULL,"-f",file,PETSC_MAX_PATH_LEN-1,&flg);CHKERRQ(ierr);
if (!flg) SETERRQ(1,"Must indicate binary file with the -f option");
ierr =
PetscViewerBinaryOpen(PETSC_COMM_WORLD,file,FILE_MODE_READ,&fd);CHKERRQ(ierr);
ierr = MatLoad(fd,MATMPIAIJ,&U);CHKERRQ(ierr);
ierr = PetscViewerDestroy(fd);CHKERRQ(ierr);
ierr = MatGetSize(U,&m,&n);CHKERRQ(ierr);
PetscPrintf(PETSC_COMM_WORLD,"# matrix Rows=%i and # matrix Columns=%i
\n\n\n",m,n);
PetscPrintf(PETSC_COMM_WORLD,"The matrix is \n\n");
// ierr=MatView(U,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
ierr=VecCreateMPI(PETSC_COMM_WORLD,PETSC_DECIDE,bsize,&b);CHKERRQ(ierr);
ierr = VecSet(b,0);CHKERRQ(ierr);
/* Reading in the RHS vector. */
fp=fopen("b1.mat","r");
if (fp==NULL)
{
fprintf(stderr, "Cannot open file");
exit(1);
}
for (i = 0; i < bsize; i++)
{
if (fscanf(fp,"%lf", &rhs[i]) != 1)
{
fprintf(stderr, "Failed to read rhs vector[%d]\n", i);
exit(1);
}
}
index=0;
/*Putting b into final form */
for (i=0; i<bsize; ++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 rhs vector. */
ierr=PetscPrintf(PETSC_COMM_WORLD,"Vector b:\n");CHKERRQ(ierr);
// ierr=VecView(b,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
/*
========================================================================================================================
*/
/* Creating the solution vector x and the new_rhs vector of size n.--->>
*/
ierr = VecCreateMPI(PETSC_COMM_WORLD,PETSC_DECIDE,n,&x);CHKERRQ(ierr);
ierr = VecCreateMPI(PETSC_COMM_WORLD,PETSC_DECIDE,n,&new_rhs);CHKERRQ(ierr);
/* Creating the U_transpose matrix and actually seting the values in it */
ierr=MatTranspose(U, MAT_INITIAL_MATRIX,&U_Transpose);CHKERRQ(ierr);
ierr=PetscPrintf(PETSC_COMM_WORLD,"\n \n U Transpose:\n\n");CHKERRQ(ierr);
// ierr=MatView(U_Transpose,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
/* Creating the matrix Prod and setting it to U_transpose*U */
ierr=MatMatMult(U_Transpose,U,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&Prod);CHKERRQ(ierr);
ierr=PetscPrintf(PETSC_COMM_WORLD,"\n\n Product Ut *U:\n");CHKERRQ(ierr);
// ierr=MatView(Prod,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
/* Creating the matrix-vec product */
ierr=MatMult(U_Transpose,b,new_rhs);CHKERRQ(ierr);
ierr=PetscPrintf(PETSC_COMM_WORLD,"\n\n new rhs:\n");CHKERRQ(ierr);
// ierr=VecView(new_rhs,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
//------------------------------------------------------
ierr = KSPCreate(PETSC_COMM_WORLD,&ksp);CHKERRQ(ierr);
/*
Set operators. Here the matrix that defines the linear system
also serves as the preconditioning matrix.
*/
ierr = KSPSetOperators(ksp,Prod,Prod,SAME_PRECONDITIONER);CHKERRQ(ierr);
ierr = KSPGetPC(ksp,&pc);CHKERRQ(ierr);
ierr = PCSetType(pc,PCJACOBI);CHKERRQ(ierr);
ierr =
KSPSetTolerances(ksp,1.e-6,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);CHKERRQ(ierr);
/*
Set runtime options, e.g.,
-ksp_type <type> -pc_type <type> -ksp_monitor -ksp_rtol <rtol>
These options will override those specified above as long as
KSPSetFromOptions() is called _after_ any other customization
routines.
*/
ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr);
ierr = KSPSolve(ksp,new_rhs,x);CHKERRQ(ierr);
ierr = KSPView(ksp,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
/*
-----------------------------------------------------------------------------------------
*/
ierr=VecView(x,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
/*
---------------------------------------------------------------------------------------
*/
ierr = MatDestroy(U);CHKERRQ(ierr);
ierr = MatDestroy(U_Transpose);CHKERRQ(ierr);
ierr = MatDestroy(Prod);CHKERRQ(ierr);
ierr = VecDestroy(b);CHKERRQ(ierr);
ierr= VecDestroy(x);CHKERRQ(ierr);
ierr= VecDestroy(new_rhs);CHKERRQ(ierr);
ierr = PetscFinalize();CHKERRQ(ierr);
return 0;
}
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20110123/ed28e0af/attachment.htm>
More information about the petsc-users
mailing list