Actual source code: test4f.F90

slepc-3.23.0 2025-03-29
Report Typos and Errors
  1: !
  2: !  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  3: !  SLEPc - Scalable Library for Eigenvalue Problem Computations
  4: !  Copyright (c) 2002-, 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:       SVDProblemType     ptype
 47:       PetscMPIInt        rank
 48:       PetscErrorCode     ierr
 49:       PetscBool          flg, tmode
 50:       PetscViewerAndFormat vf

 52: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 53: !     Beginning of program
 54: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

 56:       PetscCallA(SlepcInitialize(PETSC_NULL_CHARACTER,ierr))
 57:       PetscCallMPIA(MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr))
 58:       m = 20
 59:       PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-m',m,flg,ierr))
 60:       PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-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:       PetscCallA(MatCreate(PETSC_COMM_WORLD,A,ierr))
 73:       PetscCallA(MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,m,n,ierr))
 74:       PetscCallA(MatSetFromOptions(A,ierr))

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

 90:       PetscCallA(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY,ierr))
 91:       PetscCallA(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY,ierr))

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

 97:       PetscCallA(SVDCreate(PETSC_COMM_WORLD,svd,ierr))
 98:       PetscCallA(SVDSetOperators(svd,A,PETSC_NULL_MAT,ierr))

100: !     ** test some interface functions
101:       PetscCallA(SVDGetOperators(svd,B,PETSC_NULL_MAT,ierr))
102:       PetscCallA(MatView(B,PETSC_VIEWER_STDOUT_WORLD,ierr))
103:       PetscCallA(SVDSetConvergenceTest(svd,SVD_CONV_ABS,ierr))
104:       PetscCallA(SVDSetStoppingTest(svd,SVD_STOP_BASIC,ierr))

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

145:       PetscCallA(PetscViewerAndFormatCreate(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_DEFAULT,vf,ierr))
146:       PetscCallA(SVDMonitorSet(svd,SVDMONITORFIRST,vf,PetscViewerAndFormatDestroy,ierr))
147:       PetscCallA(SVDMonitorConvergedCreate(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_DEFAULT,PETSC_NULL_VEC,vf,ierr))
148:       PetscCallA(SVDMonitorSet(svd,SVDMONITORCONVERGED,vf,SVDMonitorConvergedDestroy,ierr))
149:       PetscCallA(SVDMonitorCancel(svd,ierr))

151: !     ** call the solver
152:       PetscCallA(SVDSetFromOptions(svd,ierr))
153:       PetscCallA(SVDSolve(svd,ierr))
154:       PetscCallA(SVDGetConvergedReason(svd,reason,ierr))
155:       if (rank .eq. 0) then
156:         write(*,150) reason
157:       endif
158:  150  format (' Converged reason:',I2)
159:       PetscCallA(SVDGetIterationNumber(svd,its,ierr))
160: !     if (rank .eq. 0) then
161: !       write(*,160) its
162: !     endif
163: !160  format (' Number of iterations of the method:',I4)

165: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
166: !     Display solution and clean up
167: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

169:       PetscCallA(SVDErrorView(svd,SVD_ERROR_RELATIVE,PETSC_NULL_VIEWER,ierr))
170:       PetscCallA(SVDDestroy(svd,ierr))
171:       PetscCallA(MatDestroy(A,ierr))

173:       PetscCallA(SlepcFinalize(ierr))
174:       end

176: !/*TEST
177: !
178: !   test:
179: !      suffix: 1
180: !      args: -svd_type {{lanczos trlanczos cross cyclic randomized}}
181: !      filter: sed -e 's/2.99255/2.99254/'
182: !
183: !TEST*/