Actual source code: test4f.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> ./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: #include <slepc/finclude/slepcsvd.h>
     29: program test4f
     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 == 0) then
     64:     write (*, '(/a,i3,a,i3,a)') 'Bidiagonal matrix, m =', m, ', n=', n, ' (Fortran)'
     65:   end if

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

     71:   PetscCallA(MatCreate(PETSC_COMM_WORLD, A, ierr))
     72:   PetscCallA(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, m, n, ierr))
     73:   PetscCallA(MatSetFromOptions(A, ierr))

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

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

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

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

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

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

    139:   PetscCallA(PetscViewerAndFormatCreate(PETSC_VIEWER_STDOUT_WORLD, PETSC_VIEWER_DEFAULT, vf, ierr))
    140:   PetscCallA(SVDMonitorSet(svd, SVDMONITORFIRST, vf, PetscViewerAndFormatDestroy, ierr))
    141:   PetscCallA(SVDMonitorConvergedCreate(PETSC_VIEWER_STDOUT_WORLD, PETSC_VIEWER_DEFAULT, PETSC_NULL, vf, ierr))
    142:   PetscCallA(SVDMonitorSet(svd, SVDMONITORCONVERGED, vf, PetscViewerAndFormatDestroy, ierr))
    143:   PetscCallA(SVDMonitorCancel(svd, ierr))

    145: ! ** call the solver
    146:   PetscCallA(SVDSetFromOptions(svd, ierr))
    147:   PetscCallA(SVDSolve(svd, ierr))
    148:   PetscCallA(SVDGetConvergedReason(svd, reason, ierr))
    149:   if (rank == 0) then
    150:     write (*, '(a,i2)') ' Converged reason:', reason
    151:   end if
    152:   PetscCallA(SVDGetIterationNumber(svd, its, ierr))
    153: ! if (rank==0) then
    154: !   write(*,'(a,i4)') ' Number of iterations of the method:', its
    155: ! end if

    157: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    158: ! Display solution and clean up
    159: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

    161:   PetscCallA(SVDErrorView(svd, SVD_ERROR_RELATIVE, PETSC_NULL_VIEWER, ierr))
    162:   PetscCallA(SVDDestroy(svd, ierr))
    163:   PetscCallA(MatDestroy(A, ierr))

    165:   PetscCallA(SlepcFinalize(ierr))
    166: end program test4f

    168: !/*TEST
    169: !
    170: !   test:
    171: !      suffix: 1
    172: !      args: -svd_type {{lanczos trlanczos cross cyclic randomized}}
    173: !      filter: sed -e 's/2.99255/2.99254/'
    174: !
    175: !TEST*/