Hi,<br><br>I have been trying to solve an overdetermined system of linear equations on PETSc using KSPLSQR. <br><br>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.)<br>
<br>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.<br><br>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. <br>
<br>Is there a way to avoid the explicit and expensive computation of U(transpose)U and use KSPLSQR? The details on this page <a href="http://www.mcs.anl.gov/petsc/petsc-2/snapshots/petsc-current/docs/manualpages/KSP/KSPLSQR.html">http://www.mcs.anl.gov/petsc/petsc-2/snapshots/petsc-current/docs/manualpages/KSP/KSPLSQR.html</a> did not really help. <br>
<br>I have also experimented with the suggesstions in the thread
<a href="http://lists.mcs.anl.gov/pipermail/petsc-users/2010-August/006784.html">http://lists.mcs.anl.gov/pipermail/petsc-users/2010-August/006784.html</a>
but I kept getting errors saying there is a dimension mismatch <br><br>I use KSPSetOperators as ierr = KSPSetOperators(ksp,Prod,Prod,SAME_PRECONDITIONER);CHKERRQ(ierr); where Prod=U(Transpose)U calculated in the previous statements.<br>
<br>
Anyway, I am providing my current working code underneath just for reference. U and b are provided through files.<br><br>I am guessing the answer lies in the way one uses KSPSetOperators but I am not sure. <br><br><br>Regards,<br>
<br>Gaurish Telang<br>%--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------<br>
<br><br><br>static char help[] = " Doing least squares \n\<br> -f <input_file> : file to load \n\n";<br><br>/* <br> Include "petscmat.h" so that we can use matrices.<br> automatically includes:<br>
petscsys.h - base PETSc routines petscvec.h - vectors<br> petscmat.h - matrices<br> petscis.h - index sets petscviewer.h - viewers <br>*/<br>#include "petscmat.h"<br>
#include "petscvec.h"<br>#include "petscksp.h" /* For the iterative solvers */<br>#include <stdlib.h><br>#include <stdio.h><br><br>#undef __FUNCT__<br>#define __FUNCT__ "main"<br>
int main(int argc,char **args)<br>{<br> Mat U,U_Transpose,Prod; /* matrix */<br> PetscViewer fd; /* viewer */<br> char file[PETSC_MAX_PATH_LEN]; /* input file name */<br>
PetscErrorCode ierr;<br> PetscTruth flg;<br> Vec b,x,new_rhs;<br> PetscInt i,n,m,bsize=2683,index; /* bsize indicates the text file size of bx.mat */<br> PetscScalar *xx;<br>
PetscScalar rhs[bsize];<br> FILE *fp;<br><br> KSP ksp;<br> PC pc;<br> PetscMPIInt size;<br><br>PetscInt num_iters;<br>PetscReal rnorm;<br>KSPConvergedReason reason;<br><br> PetscInitialize(&argc,&args,(char *)0,help);<br>
<br> ierr = PetscOptionsGetString(PETSC_NULL,"-f",file,PETSC_MAX_PATH_LEN-1,&flg);CHKERRQ(ierr);<br> if (!flg) SETERRQ(1,"Must indicate binary file with the -f option");<br> <br> ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,file,FILE_MODE_READ,&fd);CHKERRQ(ierr);<br>
ierr = MatLoad(fd,MATMPIAIJ,&U);CHKERRQ(ierr);<br> ierr = PetscViewerDestroy(fd);CHKERRQ(ierr);<br><br> ierr = MatGetSize(U,&m,&n);CHKERRQ(ierr);<br><br> PetscPrintf(PETSC_COMM_WORLD,"# matrix Rows=%i and # matrix Columns=%i \n\n\n",m,n);<br>
<br> PetscPrintf(PETSC_COMM_WORLD,"The matrix is \n\n");<br> // ierr=MatView(U,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); <br><br> ierr=VecCreateMPI(PETSC_COMM_WORLD,PETSC_DECIDE,bsize,&b);CHKERRQ(ierr);<br>
<br> ierr = VecSet(b,0);CHKERRQ(ierr);<br><br> /* Reading in the RHS vector. */<br>fp=fopen("b1.mat","r");<br>if (fp==NULL)<br>{<br> fprintf(stderr, "Cannot open file");<br> exit(1);<br>
}<br><br> for (i = 0; i < bsize; i++)<br> {<br> <br> if (fscanf(fp,"%lf", &rhs[i]) != 1)<br> {<br> fprintf(stderr, "Failed to read rhs vector[%d]\n", i);<br>
exit(1);<br> }<br> }<br> index=0;<br>/*Putting b into final form */<br>for (i=0; i<bsize; ++i)<br> {<br> ierr= VecSetValues(b,1,&index,&rhs[i],INSERT_VALUES);CHKERRQ(ierr); /* One insertion per step. Thats what the 1 in second argument stands for */<br>
index=index+1;<br> } /* The third and fourth arguments are addresses. The fifth argument is IORA */<br> /* Assembling the vector. */<br>
ierr= VecAssemblyBegin(b);CHKERRQ(ierr);<br> ierr=VecAssemblyEnd(b);CHKERRQ(ierr);<br> /* Viewing the rhs vector. */<br> ierr=PetscPrintf(PETSC_COMM_WORLD,"Vector b:\n");CHKERRQ(ierr);<br> // ierr=VecView(b,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); <br>
<br>/* ======================================================================================================================== */<br> <br><br> /* Creating the solution vector x and the new_rhs vector of size n.--->> */<br>
ierr = VecCreateMPI(PETSC_COMM_WORLD,PETSC_DECIDE,n,&x);CHKERRQ(ierr);<br>ierr = VecCreateMPI(PETSC_COMM_WORLD,PETSC_DECIDE,n,&new_rhs);CHKERRQ(ierr);<br> /* Creating the U_transpose matrix and actually seting the values in it */<br>
<br> ierr=MatTranspose(U, MAT_INITIAL_MATRIX,&U_Transpose);CHKERRQ(ierr);<br> ierr=PetscPrintf(PETSC_COMM_WORLD,"\n \n U Transpose:\n\n");CHKERRQ(ierr);<br> // ierr=MatView(U_Transpose,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); <br>
<br> /* Creating the matrix Prod and setting it to U_transpose*U */ <br><br> ierr=MatMatMult(U_Transpose,U,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&Prod);CHKERRQ(ierr);<br> ierr=PetscPrintf(PETSC_COMM_WORLD,"\n\n Product Ut *U:\n");CHKERRQ(ierr);<br>
// ierr=MatView(Prod,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); <br><br> /* Creating the matrix-vec product */<br> ierr=MatMult(U_Transpose,b,new_rhs);CHKERRQ(ierr);<br> ierr=PetscPrintf(PETSC_COMM_WORLD,"\n\n new rhs:\n");CHKERRQ(ierr);<br>
// ierr=VecView(new_rhs,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); <br><br> //------------------------------------------------------<br> ierr = KSPCreate(PETSC_COMM_WORLD,&ksp);CHKERRQ(ierr);<br><br> /* <br> Set operators. Here the matrix that defines the linear system<br>
also serves as the preconditioning matrix.<br> */<br> ierr = KSPSetOperators(ksp,Prod,Prod,SAME_PRECONDITIONER);CHKERRQ(ierr);<br><br> ierr = KSPGetPC(ksp,&pc);CHKERRQ(ierr);<br> ierr = PCSetType(pc,PCJACOBI);CHKERRQ(ierr);<br>
ierr = KSPSetTolerances(ksp,1.e-6,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);CHKERRQ(ierr);<br><br> /* <br> Set runtime options, e.g.,<br> -ksp_type <type> -pc_type <type> -ksp_monitor -ksp_rtol <rtol><br>
These options will override those specified above as long as<br> KSPSetFromOptions() is called _after_ any other customization<br> routines.<br> */<br> ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr);<br> <br> ierr = KSPSolve(ksp,new_rhs,x);CHKERRQ(ierr); <br>
<br><br> ierr = KSPView(ksp,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);<br><br> /* ----------------------------------------------------------------------------------------- */<br><br> ierr=VecView(x,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); <br>
<br> /* --------------------------------------------------------------------------------------- */<br> ierr = MatDestroy(U);CHKERRQ(ierr);<br> ierr = MatDestroy(U_Transpose);CHKERRQ(ierr);<br> ierr = MatDestroy(Prod);CHKERRQ(ierr); <br>
<br> ierr = VecDestroy(b);CHKERRQ(ierr);<br> ierr= VecDestroy(x);CHKERRQ(ierr); <br> ierr= VecDestroy(new_rhs);CHKERRQ(ierr); <br><br>ierr = PetscFinalize();CHKERRQ(ierr);<br> return 0;<br>}<br><br><br><br><br>