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>
