Actual source code: ex15f.F90

      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> ./ex15f [-help] [-n <n>] [-mu <mu>] [all SLEPc options]
     11: !
     12: !  Description: Singular value decomposition of the Lauchli matrix.
     13: !
     14: !  The command line options are:
     15: !    -n <n>, where <n> = matrix dimension.
     16: !    -mu <mu>, where <mu> = subdiagonal value.
     17: !
     18: ! ----------------------------------------------------------------------
     19: !
     20: #include <slepc/finclude/slepcsvd.h>
     21: program ex15f
     22:   use slepcsvd
     23:   implicit none

     25: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     26: ! Declarations
     27: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

     29:   Mat            :: A    ! operator matrix
     30:   SVD            :: svd  ! singular value solver context
     31:   SVDType        :: tname
     32:   PetscReal      :: tol, error, sigma, mu
     33:   PetscInt       :: n, i, j, Istart, Iend
     34:   PetscInt       :: nsv, maxit, its, nconv
     35:   PetscMPIInt    :: rank
     36:   PetscErrorCode :: ierr
     37:   PetscBool      :: flg
     38:   PetscScalar    :: alpha
     39:   PetscScalar, parameter :: one = 1.0

     41: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     42: ! Beginning of program
     43: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

     45:   PetscCallA(SlepcInitialize(PETSC_NULL_CHARACTER, ierr))
     46:   PetscCallMPIA(MPI_Comm_rank(PETSC_COMM_WORLD, rank, ierr))
     47:   n = 100
     48:   PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-n', n, flg, ierr))
     49:   mu = PETSC_SQRT_MACHINE_EPSILON
     50:   PetscCallA(PetscOptionsGetReal(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-mu', mu, flg, ierr))

     52:   if (rank == 0) then
     53:     write (*, '(/a,i3,a,e12.4,a)') 'Lauchli SVD, n =', n, ', mu=', mu, ' (Fortran)'
     54:   end if

     56: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     57: ! Build the Lauchli matrix
     58: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

     60:   PetscCallA(MatCreate(PETSC_COMM_WORLD, A, ierr))
     61:   PetscCallA(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, n + 1, n, ierr))
     62:   PetscCallA(MatSetFromOptions(A, ierr))

     64:   PetscCallA(MatGetOwnershipRange(A, Istart, Iend, ierr))
     65:   do i = Istart, Iend - 1
     66:     if (i == 0) then
     67:       do j = 0, n - 1
     68:         PetscCallA(MatSetValue(A, i, j, one, INSERT_VALUES, ierr))
     69:       end do
     70:     else
     71:       alpha = mu
     72:       PetscCallA(MatSetValue(A, i, i - 1, alpha, INSERT_VALUES, ierr))
     73:     end if
     74:   end do

     76:   PetscCallA(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY, ierr))
     77:   PetscCallA(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY, ierr))

     79: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     80: ! Create the singular value solver and display info
     81: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

     83: ! ** Create singular value solver context
     84:   PetscCallA(SVDCreate(PETSC_COMM_WORLD, svd, ierr))

     86: ! ** Set operators and problem type
     87:   PetscCallA(SVDSetOperators(svd, A, PETSC_NULL_MAT, ierr))
     88:   PetscCallA(SVDSetProblemType(svd, SVD_STANDARD, ierr))

     90: ! ** Use thick-restart Lanczos as default solver
     91:   PetscCallA(SVDSetType(svd, SVDTRLANCZOS, ierr))

     93: ! ** Set solver parameters at runtime
     94:   PetscCallA(SVDSetFromOptions(svd, ierr))

     96: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     97: ! Solve the singular value system
     98: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

    100:   PetscCallA(SVDSolve(svd, ierr))
    101:   PetscCallA(SVDGetIterationNumber(svd, its, ierr))
    102:   if (rank == 0) then
    103:     write (*, '(/a,i4)') ' Number of iterations of the method:', its
    104:   end if

    106: ! ** Optional: Get some information from the solver and display it
    107:   PetscCallA(SVDGetType(svd, tname, ierr))
    108:   if (rank == 0) then
    109:     write (*, '(a,a)') ' Solution method: ', tname
    110:   end if
    111:   PetscCallA(SVDGetDimensions(svd, nsv, PETSC_NULL_INTEGER, PETSC_NULL_INTEGER, ierr))
    112:   if (rank == 0) then
    113:     write (*, '(a,i2)') ' Number of requested singular values:', nsv
    114:   end if
    115:   PetscCallA(SVDGetTolerances(svd, tol, maxit, ierr))
    116:   if (rank == 0) then
    117:     write (*, '(a,1pe11.4,a,i4)') ' Stopping condition: tol=', tol, ', maxit=', maxit
    118:   end if

    120: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    121: ! Display solution and clean up
    122: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

    124: ! ** Get number of converged singular triplets
    125:   PetscCallA(SVDGetConverged(svd, nconv, ierr))
    126:   if (rank == 0) then
    127:     write (*, '(a,i2/)') ' Number of converged approximate singular triplets:', nconv
    128:   end if

    130: ! ** Display singular values and relative errors
    131:   if (nconv > 0) then
    132:     if (rank == 0) then
    133:       write (*, *) '       sigma          relative error'
    134:       write (*, *) ' ----------------- ------------------'
    135:     end if
    136:     do i = 0, nconv - 1
    137: !     ** Get i-th singular value
    138:       PetscCallA(SVDGetSingularTriplet(svd, i, sigma, PETSC_NULL_VEC, PETSC_NULL_VEC, ierr))

    140: !     ** Compute the relative error for each singular triplet
    141:       PetscCallA(SVDComputeError(svd, i, SVD_ERROR_RELATIVE, error, ierr))
    142:       if (rank == 0) then
    143:         write (*, '(1p,a,e12.4,a,e12.4)') '   ', sigma, '       ', error
    144:       end if

    146:     end do
    147:     if (rank == 0) then
    148:       write (*, *)
    149:     end if
    150:   end if

    152: ! ** Free work space
    153:   PetscCallA(SVDDestroy(svd, ierr))
    154:   PetscCallA(MatDestroy(A, ierr))

    156:   PetscCallA(SlepcFinalize(ierr))
    157: end program ex15f

    159: !/*TEST
    160: !
    161: !   test:
    162: !      suffix: 1
    163: !      filter: sed -e "s/[0-9]\.[0-9]*E[+-]\([0-9]*\)/removed/g"
    164: !
    165: !TEST*/