Actual source code: test4f.F

slepc-3.11.2 2019-07-30
Report Typos and Errors
  1: !
  2: !  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  3: !  SLEPc - Scalable Library for Eigenvalue Problem Computations
  4: !  Copyright (c) 2002-2019, Universitat Politecnica de Valencia, Spain
  5: !
  6: !  This file is part of SLEPc.
  7: !  SLEPc is distributed under a 2-clause BSD license (see LICENSE).
  8: !  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  9: !
 10: !  Program usage: mpiexec -n <np> ./test4f [-help] [-n <n>] [-m <m>] [all SLEPc options]
 11: !
 12: !  Description: Singular value decomposition of a bidiagonal matrix.
 13: !
 14: !               |  1  2                     |
 15: !               |     1  2                  |
 16: !               |        1  2               |
 17: !           A = |          .  .             |
 18: !               |             .  .          |
 19: !               |                1  2       |
 20: !               |                   1  2    |
 21: !
 22: !  The command line options are:
 23: !    -m <m>, where <m> = matrix rows.
 24: !    -n <n>, where <n> = matrix columns (defaults to m+2).
 25: !
 26: ! ----------------------------------------------------------------------
 27: !
 28:       program main
 29: #include <slepc/finclude/slepcsvd.h>
 30:       use slepcsvd
 31:       implicit none

 33: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 34: !     Declarations
 35: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 36: !
 37:       Mat                A, B
 38:       SVD                svd
 39:       SVDConv            conv;
 40:       SVDStop            stp;
 41:       SVDWhich           which;
 42:       SVDConvergedReason reason;
 43:       PetscInt           m, n, i, Istart
 44:       PetscInt           col(2), its, Iend
 45:       PetscScalar        val(2)
 46:       PetscMPIInt        rank
 47:       PetscErrorCode     ierr
 48:       PetscBool          flg, tmode

 50: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 51: !     Beginning of program
 52: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

 54:       call SlepcInitialize(PETSC_NULL_CHARACTER,ierr)
 55:       call MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr)
 56:       m = 20
 57:       call PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,  &
 58:      &                        '-m',m,flg,ierr)
 59:       call PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,  &
 60:      &                        '-n',n,flg,ierr)
 61:       if (.not. flg) n = m+2

 63:       if (rank .eq. 0) then
 64:         write(*,100) m, n
 65:       endif
 66:  100  format (/'Bidiagonal matrix, m =',I3,', n=',I3,' (Fortran)')

 68: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 69: !     Build the Lauchli matrix
 70: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

 72:       call MatCreate(PETSC_COMM_WORLD,A,ierr)
 73:       call MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,m,n,ierr)
 74:       call MatSetFromOptions(A,ierr)
 75:       call MatSetUp(A,ierr)

 77:       call MatGetOwnershipRange(A,Istart,Iend,ierr)
 78:       val(1) = 1.0
 79:       val(2) = 2.0
 80:       do i=Istart,Iend-1
 81:         col(1) = i
 82:         col(2) = i+1
 83:         if (i .le. n-1) then
 84:           call MatSetValue(A,i,col(1),val(1),INSERT_VALUES,ierr)
 85:         end if
 86:         if (i .lt. n-1) then
 87:           call MatSetValue(A,i,col(2),val(2),INSERT_VALUES,ierr)
 88:         end if
 89:       enddo

 91:       call MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY,ierr)
 92:       call MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY,ierr)

 94: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 95: !     Compute singular values
 96: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

 98:       call SVDCreate(PETSC_COMM_WORLD,svd,ierr)
 99:       call SVDSetOperator(svd,A,ierr)

101: !     ** test some interface functions
102:       call SVDGetOperator(svd,B,ierr)
103:       call MatView(B,PETSC_VIEWER_STDOUT_WORLD,ierr)
104:       call SVDSetConvergenceTest(svd,SVD_CONV_ABS,ierr)
105:       call SVDSetStoppingTest(svd,SVD_STOP_BASIC,ierr)
106:       call SVDSetFromOptions(svd,ierr)

108: !     ** query properties and print them
109:       call SVDGetImplicitTranspose(svd,tmode,ierr)
110:       if (rank .eq. 0) then
111:         if (tmode) then
112:           write(*,110) 'implicit'
113:         else
114:           write(*,110) 'explicit'
115:         endif
116:       endif
117:  110  format (/' Transpose mode is',A9)
118:       call SVDGetConvergenceTest(svd,conv,ierr)
119:       if (rank .eq. 0) then
120:         write(*,120) conv
121:       endif
122:  120  format (' Convergence test is',I2)
123:       call SVDGetStoppingTest(svd,stp,ierr)
124:       if (rank .eq. 0) then
125:         write(*,130) stp
126:       endif
127:  130  format (' Stopping test is',I2)
128:       call SVDGetWhichSingularTriplets(svd,which,ierr)
129:       if (rank .eq. 0) then
130:         if (which .eq. SVD_LARGEST) then
131:           write(*,140) 'largest'
132:         else
133:           write(*,140) 'smallest'
134:         endif
135:       endif
136:  140  format (' Which =',A9)

138: !     ** call the solver
139:       call SVDSolve(svd,ierr)
140:       call SVDGetConvergedReason(svd,reason,ierr)
141:       if (rank .eq. 0) then
142:         write(*,150) reason
143:       endif
144:  150  format (' Converged reason:',I2)
145:       call SVDGetIterationNumber(svd,its,ierr)
146: !     if (rank .eq. 0) then
147: !       write(*,160) its
148: !     endif
149: !160  format (' Number of iterations of the method:',I4)

151: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
152: !     Display solution and clean up
153: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

155:       call SVDErrorView(svd,SVD_ERROR_RELATIVE,PETSC_NULL_VIEWER,ierr)
156:       call SVDDestroy(svd,ierr)
157:       call MatDestroy(A,ierr)

159:       call SlepcFinalize(ierr)
160:       end

162: !/*TEST
163: !
164: !   test:
165: !      suffix: 1
166: !
167: !TEST*/