CoCalc -- EigenSolve.F90

GitHub Repository: ElmerCSC/elmerfem
Path: blob/devel/fem/src/EigenSolve.F90
⁵²²⁶ views
1
!/*****************************************************************************/
2
! *
3
! *  Elmer, A Finite Element Software for Multiphysical Problems
4
! *
5
! *  Copyright 1st April 1995 - , CSC - IT Center for Science Ltd., Finland
6
! * 
7
! *  This library is free software; you can redistribute it and/or
8
! *  modify it under the terms of the GNU Lesser General Public
9
! *  License as published by the Free Software Foundation; either
10
! *  version 2.1 of the License, or (at your option) any later version.
11
! *
12
! *  This library is distributed in the hope that it will be useful,
13
! *  but WITHOUT ANY WARRANTY; without even the implied warranty of
14
! *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
15
! *  Lesser General Public License for more details.
16
! * 
17
! *  You should have received a copy of the GNU Lesser General Public
18
! *  License along with this library (in file ../LGPL-2.1); if not, write 
19
! *  to the Free Software Foundation, Inc., 51 Franklin Street, 
20
! *  Fifth Floor, Boston, MA  02110-1301  USA
21
! *
22
! *****************************************************************************/
23
!
24
!/******************************************************************************
25
! *
26
! *  This module provides various solution techniques for eigenvalue problems. 
27
! *  The module is based on the ARPACK example driver dndrv3 that has been modified
28
! *  to fit the needs of ElmerSolver.
29
! *
30
! *
31
! *  FILE: ndrv3.F   SID: 2.4   DATE OF SID: 4/22/96   RELEASE: 2
32
! *
33
! *  Example Authors: Richard Lehoucq, Danny Sorensen and Chao Yang
34
! *                   Dept. of Computational & Applied Mathematics
35
! *                   Rice University
36
! *                   Houston, Texas
37
! *
38
! ******************************************************************************
39
! *
40
! *  Authors: Juha Ruokolainen
41
! *  Email:   [email protected]
42
! *  Web:     http://www.csc.fi/elmer
43
! *  Address: CSC - IT Center for Science Ltd.
44
! *           Keilaranta 14
45
! *           02101 Espoo, Finland 
46
! *
47
! *  Original Date: 21 Oct 2000
48
! *
49
! *****************************************************************************/
50

51

52
!----------------------------------------------------------------------------
53
!> Module containing ARPACK routines for the solution of Eigenvalue problems.
54
!----------------------------------------------------------------------------
55
!> \ingroup ElmerLib 
56
!> \{
57

58
MODULE EigenSolve
59

60
   IMPLICIT NONE
61

62
CONTAINS
63

64
!------------------------------------------------------------------------------
65
!> Solution of Eigen value problems using ARPACK library. 
66
!------------------------------------------------------------------------------
67
     SUBROUTINE ArpackEigenSolve( Solver,Matrix,N,NEIG,EigValues,EigVectors )
68
!------------------------------------------------------------------------------
69
      USE Multigrid
70

71
      IMPLICIT NONE
72

73
      TYPE(Matrix_t), POINTER :: Matrix, A
74
      TYPE(Solver_t), TARGET :: Solver
75
      INTEGER :: N, NEIG
76
      INTEGER, ALLOCATABLE :: dperm(:)
77
      COMPLEX(KIND=dp) :: EigValues(:), EigVectors(:,:)
78

79
#ifdef USE_ARPACK
80
!
81
!     %--------------%
82
!     | Local Arrays |
83
!     %--------------%
84
!
85
      REAL(KIND=dp), TARGET, ALLOCATABLE :: WORKD(:), RESID(:),bb(:),xx(:)
86
      REAL(KIND=dp), POINTER CONTIG :: x(:), b(:)
87
      INTEGER :: IPARAM(11), IPNTR(14)
88
      INTEGER, ALLOCATABLE :: Perm(:)
89
      LOGICAL, ALLOCATABLE :: Choose(:)
90
      CHARACTER(:), ALLOCATABLE :: Method
91
      REAL(KIND=dp), ALLOCATABLE :: WORKL(:), D(:,:), WORKEV(:), V(:,:)
92

93
!
94
!     %---------------%
95
!     | Local Scalars |
96
!     %---------------%
97
!
98
      CHARACTER ::     BMAT*1, Which*2, DirectMethod*100
99
      INTEGER   ::     IDO, NCV, lWORKL, kinfo, i, j, k, l, p, IERR, iter, &
100
                       NCONV, maxitr, ishfts, mode, istat, Dofs
101
      LOGICAL   ::     First, Stat, Direct = .FALSE., &
102
                       Iterative = .FALSE., NewSystem, Damped, Stability, &
103
                       NormalizeToUnity
104

105
      LOGICAL :: Factorize, FreeFactorize,FoundFactorize,FoundFreeFactorize
106
      REAL(KIND=dp) :: SigmaR, SigmaI, TOL, r
107

108
      COMPLEX(KIND=dp) :: s
109
!
110
      REAL(KIND=dp), POINTER CONTIG :: SaveValues(:), SaveRhs(:)
111
      TYPE(ValueList_t), POINTER :: Params
112

113
      CHARACTER(*), PARAMETER :: Caller = 'EigenSolve'
114

115
      
116
!     %--------------------------------------%
117
!     | Check if system is damped and if so, |
118
!     | move to other subroutine             |
119
!     %--------------------------------------%
120

121
      Params => Solver % Values
122
      Damped = ListGetLogical( Params, 'Eigen System Damped', stat )
123
      IF ( .NOT. stat ) Damped = .FALSE.
124

125
      IF ( Damped ) THEN
126
         CALL ArpackDampedEigenSolve( Solver, Matrix, 2*N, 2*NEIG, &
127
              EigValues,EigVectors )
128
         RETURN
129
      END IF
130

131
!     %----------------------------------------%
132
!     | Check if stability analysis is defined |
133
!     | and if so move to other subroutine     |
134
!     %----------------------------------------%
135

136
      Stability = ListGetLogical( Params, 'stability analysis', stat )
137
      IF ( .NOT. stat ) Stability = .FALSE.
138

139
      IF ( Stability ) THEN
140
         CALL ArpackStabEigenSolve( Solver, Matrix, N, NEIG, EigValues,EigVectors )
141
         RETURN
142
      END IF
143

144
!     %-----------------------%
145
!     | Executable Statements |
146
!     %-----------------------%
147
!
148
!     %----------------------------------------------------%
149
!     | The number N is the dimension of the matrix. A     |
150
!     | generalized eigenvalue problem is solved (BMAT =   |
151
!     | 'G'.) NEV is the number of eigenvalues to be       |
152
!     | approximated.  The user can modify NEV, NCV, WHICH |
153
!     | to solve problems of different sizes, and to get   |
154
!     | different parts of the spectrum.  However, The     |
155
!     | following conditions must be satisfied:            |
156
!     |                     N <= MAXN,                     | 
157
!     |                   NEV <= MAXNEV,                   |
158
!     |               NEV + 1 <= NCV <= MAXNCV             | 
159
!     %----------------------------------------------------%
160
!
161
      NCV = ListGetInteger( Params, 'Eigen System Lanczos Vectors', stat )
162
      IF ( .NOT. stat ) NCV = 3*NEIG + 1
163

164
      IF ( NCV <=  NEIG ) THEN
165
         CALL Fatal( Caller, & 
166
               'Number of Lanczos vectors must exceed the number of eigenvalues.' )
167
      END IF
168

169
      ALLOCATE(workd(3*n),resid(n), bb(n), xx(n), stat=istat )
170
      IF ( istat /= 0 ) CALL Fatal( Caller, 'Memory allocation error.' )
171

172
      ALLOCATE( WORKL(3*NCV**2 + 6*NCV), D(NCV,3), &
173
                WORKEV(3*NCV), V(n,NCV), CHOOSE(NCV), STAT=istat )
174
      IF ( istat /= 0 ) CALL Fatal( Caller, 'Memory allocation error.' )
175

176
!
177
!     %--------------------------------------------------%
178
!     | The work array WORKL is used in DSAUPD as        |
179
!     | workspace.  Its dimension LWORKL is set as       |
180
!     | illustrated below.  The parameter TOL determines |
181
!     | the stopping criterion.  If TOL<=0, machine      |
182
!     | precision is used.  The variable IDO is used for |
183
!     | reverse communication and is initially set to 0. |
184
!     | Setting INFO=0 indicates that a random vector is |
185
!     | generated in DSAUPD to start the Arnoldi         |
186
!     | iteration.                                       |
187
!     %--------------------------------------------------%
188
!
189
!
190

191
      TOL = ListGetConstReal( Params, 'Eigen System Convergence Tolerance', stat )
192
      IF ( .NOT. stat ) THEN
193
         TOL = 100 * ListGetConstReal( Params, 'Linear System Convergence Tolerance' )
194
      END IF
195

196
      IDO   = 0
197
      kinfo = 0
198
      lWORKL = 3*NCV**2 + 6*NCV 
199
!
200
!     %---------------------------------------------------%
201
!     | This program uses exact shifts with respect to    |
202
!     | the current Hessenberg matrix (IPARAM(1) = 1).    |
203
!     | IPARAM(3) specifies the maximum number of Arnoldi |
204
!     | iterations allowed.  Mode 2 of DSAUPD is used     |
205
!     | (IPARAM(7) = 2).  All these options may be        |
206
!     | changed by the user. For details, see the         |
207
!     | documentation in DSAUPD.                          |
208
!     %---------------------------------------------------%
209
!
210
      ishfts = 1
211
      BMAT  = 'G'
212
      IF ( Matrix % Lumped ) THEN
213
         Mode  =  2
214
         SELECT CASE( ListGetString(Params,'Eigen System Select',stat) )
215
         CASE( 'smallest magnitude' )
216
              Which = 'SM'
217
         CASE( 'largest magnitude')
218
              Which = 'LM'
219
         CASE( 'smallest real part')
220
              Which = 'SR'
221
         CASE( 'largest real part')
222
              Which = 'LR'
223
         CASE( 'smallest imag part' )
224
              Which = 'SI'
225
         CASE( 'largest imag part' )
226
              Which = 'LI'
227
         CASE DEFAULT
228
              Which = 'SM'
229
         END SELECT
230
      ELSE
231
         Mode  = 3
232
         SELECT CASE( ListGetString(Params,'Eigen System Select',stat) )
233
         CASE( 'smallest magnitude' )
234
              Which = 'LM'
235
         CASE( 'largest magnitude')
236
              Which = 'SM'
237
         CASE( 'smallest real part')
238
              Which = 'LR'
239
         CASE( 'largest real part')
240
              Which = 'SR'
241
         CASE( 'smallest imag part' )
242
              Which = 'LI'
243
         CASE( 'largest imag part' )
244
              Which = 'SI'
245
         CASE DEFAULT
246
              Which = 'LM'
247
         END SELECT
248
      END IF
249

250
      Maxitr = ListGetInteger( Params, 'Eigen System Max Iterations', stat )
251
      IF ( .NOT. stat ) Maxitr = 300
252
!
253
      IPARAM = 0
254
      IPARAM(1) = ishfts
255
      IPARAM(3) = maxitr 
256
      IPARAM(7) = mode
257

258
      SigmaR = 0.0d0
259
      SigmaI = 0.0d0
260
      V = 0.0d0
261

262
!     Compute LU-factors for (A-\sigma M) (if consistent mass matrix)
263
!
264
      Factorize = ListGetLogical( Params, &
265
            'Linear System Refactorize', FoundFactorize )
266
      CALL ListAddLogical( Params, 'Linear System Refactorize',.TRUE. )
267

268
      FreeFactorize = ListGetLogical( Params, &
269
                'Linear System Free Fctorization', FoundFreeFactorize )
270
      CALL ListAddLogical( Params,  &
271
                     'Linear System Free Factorization',.FALSE. )
272

273
      IF ( .NOT. Matrix % Lumped ) THEN
274
        SigmaR = ListGetConstReal( Params,'Eigen System Shift', stat )
275
        IF ( SigmaR /= 0.0d0 ) THEN
276
          Matrix % Values = Matrix % Values - SigmaR * Matrix % MassValues
277
        END IF
278

279
        Method = ListGetString( Params,'Linear System Solver', stat )         
280
        IF ( Method == 'direct' ) THEN
281
          DirectMethod = ListGetString( Params, &
282
              'Linear System Direct Method', stat )
283
          
284
          SELECT CASE( DirectMethod )
285
          CASE('umfpack', 'big umfpack', 'mumps', 'superlu', 'pardiso', 'cholmod')
286
          CASE DEFAULT
287
            Stat = CRS_ILUT(Matrix, 0.0d0)
288
          END SELECT
289
        END IF
290
      END IF
291

292
!
293
!     %-------------------------------------------%
294
!     | M A I N   L O O P (Reverse communication) |
295
!     %-------------------------------------------%
296
!
297
      iter = 1
298
      NewSystem = .TRUE.
299

300
      Iterative = ListGetString( Params, &
301
               'Linear System Solver', stat ) == 'iterative'
302

303
      stat = ListGetLogical( Params,  'No Precondition Recompute', stat  )
304
      IF ( Iterative .AND. Stat ) THEN
305
         CALL ListAddLogical( Params, 'No Precondition Recompute', .FALSE. )
306
      END IF
307

308
      DO WHILE( ido /= 99 )
309

310
!        %---------------------------------------------%
311
!        | Repeatedly call the routine DSAUPD and take | 
312
!        | actions indicated by parameter IDO until    |
313
!        | either convergence is indicated or maxitr   |
314
!        | has been exceeded.                          |
315
!        %---------------------------------------------%
316
!
317
         IF ( Matrix % Symmetric ) THEN
318
            CALL DSAUPD ( ido, BMAT, n, Which, NEIG, TOL, &
319
              RESID, NCV, V, n, IPARAM, IPNTR, WORKD, WORKL, lWORKL, kinfo )
320
         ELSE
321
            CALL DNAUPD ( ido, BMAT, n, Which, NEIG, TOL, &
322
              RESID, NCV, v, n, IPARAM, IPNTR, WORKD, WORKL, lWORKL, kinfo )
323
         END IF
324
 
325
         IF (ido == -1 .OR. ido == 1) THEN
326

327
            IF(InfoActive(20)) THEN
328
              CALL Info(Caller,'Arpack reverse communication calls: '//I2S(iter))
329
            ELSE
330
              CALL Info(Caller, '.', .TRUE.)
331
            END IF
332
            iter = iter + 1
333

334
!---------------------------------------------------------------------
335
!             Perform  y <--- OP*x = inv[M]*A*x   (lumped mass)
336
!                      ido =-1 inv(A-sigmaR*M)*M*x 
337
!                      ido = 1 inv(A-sigmaR*M)*z
338
!---------------------------------------------------------------------
339

340
            IF ( Matrix % Lumped ) THEN
341
              CALL CRS_MatrixVectorMultiply( Matrix, WORKD(IPNTR(1)), WORKD(IPNTR(2)) )
342
              DO i=0,n-1
343
                WORKD( IPNTR(1)+i ) = WORKD( IPNTR(2)+i )
344
              END DO
345
              DO i=0,n-1
346
                WORKD( IPNTR(2)+i ) = WORKD( IPNTR(1)+i ) / &
347
                    Matrix % MassValues( Matrix % Diag(i+1) )
348
              END DO
349
            ELSE              
350
             
351
              Dofs = Solver % Variable % Dofs 
352
              A => Matrix
353
              x => workd(ipntr(2):ipntr(2)+n-1)
354

355
              IF ( ido == -1 ) THEN
356
                SaveValues => A % Values
357
                A % Values => A % MassValues
358
                CALL CRS_MatrixVectorMultiply( A, WORKD(IPNTR(1)), WORKD(IPNTR(2)) )
359
                A % Values => SaveValues
360
                
361
                DO i=0,n-1
362
                  WORKD( IPNTR(1)+i ) = WORKD( IPNTR(2)+i )
363
                END DO
364
                b => workd(ipntr(1):ipntr(1)+n-1)
365
              ELSE
366
                b => workd(ipntr(3):ipntr(3)+n-1)                              
367
              END IF
368

369
              ! Some strategies (such as 'block') may depend on that these are set properly 
370
              ! to reflect the linear problem under study.            
371
              SaveRhs => A % rhs
372
              A % rhs => b
373

374
              SELECT CASE( Method ) 
375
              CASE('multigrid')
376
                CALL MultiGridSolve( A, x, b, &
377
                    DOFs, Solver, Solver % MultiGridLevel, NewSystem )
378
              CASE('iterative')
379
                CALL IterSolver( A, x, b, Solver )
380
              CASE('block')
381
                CALL BlockSolveExt( A, x, b, Solver )
382
              CASE ('direct')
383
                CALL DirectSolver( A, x, b, Solver )
384
              CASE DEFAULT
385
                CALL Fatal(Caller,'Unknown linear system method: '//TRIM(Method))
386
              END SELECT
387

388
              A % rhs => SaveRhs
389
            END IF
390

391
          ELSE IF (ido == 2) THEN
392

393
!
394
!           %-----------------------------------------%
395
!           |         Perform  y <--- M*x.            |
396
!           | Need the matrix vector multiplication   |
397
!           | routine here that takes WORKD(IPNTR(1)) |
398
!           | as the input and returns the result to  |
399
!           | WORKD(IPNTR(2)).                        |
400
!           %-----------------------------------------%
401
!
402
            IF ( Matrix % Lumped ) THEN
403
               DO i=0,n-1
404
                  WORKD( IPNTR(2)+i ) = WORKD( IPNTR(1)+i ) * &
405
                    Matrix % MassValues( Matrix % Diag(i+1) )
406
               END DO
407
            ELSE
408
               SaveValues => Matrix % Values
409
               Matrix % Values => Matrix % MassValues
410
               CALL CRS_MatrixVectorMultiply( Matrix, WORKD(IPNTR(1)), WORKD(IPNTR(2)) )
411
               Matrix % Values => SaveValues
412
            END IF
413
         END IF 
414

415
         IF ( NewSystem .AND. ido /= 2 ) THEN
416
            IF ( Iterative ) THEN
417
               CALL ListAddLogical( Params,  'No Precondition Recompute', .TRUE. )
418
            ELSE
419
               CALL ListAddLogical( Params, 'Linear System Refactorize', .FALSE. )
420
            END IF
421
            NewSystem = .FALSE.
422
         END IF
423
       END DO
424

425
      IF ( FoundFactorize ) THEN
426
        CALL ListAddLogical( Params, 'Linear System Refactorize', Factorize )
427
      ELSE
428
        CALL ListRemove( Params, 'Linear System Refactorize' )
429
      END IF
430

431
      IF ( .NOT. FoundFreeFactorize ) THEN
432
        CALL ListRemove( Params, 'Linear System Free Factorization' )
433
      ELSE
434
        CALL ListAddLogical( Params, 'Linear System Free Factorization', FreeFactorize )
435
      END IF
436

437
!     %-----------------------------------------%
438
!     | Either we have convergence, or there is |
439
!     | an error.                               |
440
!     %-----------------------------------------%
441
      IF ( kinfo /= 0 ) THEN
442
!
443
!        %--------------------------%
444
!        | Error message, check the |
445
!        | documentation in DNAUPD  |
446
!        %--------------------------%
447
!
448
         WRITE( Message, * ) 'Error with DNAUPD, info = ',kinfo
449
         CALL Fatal( Caller, Message )
450
!
451
      ELSE 
452
!
453
!        %-------------------------------------------%
454
!        | No fatal errors occurred.                 |
455
!        | Post-Process using DSEUPD.                |
456
!        |                                           |
457
!        | Computed eigenvalues may be extracted.    |  
458
!        |                                           |
459
!        | Eigenvectors may also be computed now if  |
460
!        | desired.  (indicated by rvec = .true.)    | 
461
!        %-------------------------------------------%
462
!           
463
         D = 0.0d0
464
         IF ( Matrix % Symmetric ) THEN
465
            CALL DSEUPD ( .TRUE., 'A', Choose, D, V, N, SigmaR,  &
466
               BMAT, n, Which, NEIG, TOL, RESID, NCV, V, N, &
467
               IPARAM, IPNTR, WORKD, WORKL, lWORKL, IERR )
468
         ELSE
469
            CALL DNEUPD ( .TRUE., 'A', Choose, D, D(1,2), &
470
               V, N, SigmaR, SigmaI, WORKEV, BMAT, N, &
471
               Which, NEIG, TOL, RESID, NCV, V, N, &
472
               IPARAM, IPNTR, WORKD, WORKL, lWORKL, IERR )
473
         END IF
474
 
475
!        %----------------------------------------------%
476
!        | Eigenvalues are returned in the First column |
477
!        | of the two dimensional array D and the       |
478
!        | corresponding eigenvectors are returned in   |
479
!        | the First NEV columns of the two dimensional |
480
!        | array V if requested.  Otherwise, an         |
481
!        | orthogonal basis for the invariant subspace  |
482
!        | corresponding to the eigenvalues in D is     |
483
!        | returned in V.                               |
484
!        %----------------------------------------------%
485
!
486
         IF (IERR /= 0) THEN 
487
!
488
!           %------------------------------------%
489
!           | Error condition:                   |
490
!           | Check the documentation of DNEUPD. |
491
!           %------------------------------------%
492
! 
493
            WRITE( Message, * ) ' Error with DNEUPD, info = ', IERR
494
            CALL Fatal( Caller, Message )
495
         END IF
496
!
497
!        %------------------------------------------%
498
!        | Print additional convergence information |
499
!        %------------------------------------------%
500
!
501
         IF ( kinfo == 1 ) THEN
502
            CALL Fatal( Caller, 'Maximum number of iterations reached.' )
503
         ELSE IF ( kinfo == 3 ) THEN
504
            CALL Fatal( Caller, 'No shifts could be applied during implicit Arnoldi update, try increasing NCV.' )
505
         END IF      
506
!
507
!        Sort the eigenvalues to ascending order:
508
!        ----------------------------------------
509
         ALLOCATE( Perm(NEIG) )
510
         Perm = [ (i, i=1,NEIG) ]
511
         DO i=1,NEIG
512
            EigValues(i) = CMPLX( D(i,1), D(i,2),KIND=dp )
513
         END DO
514
         CALL SortC( NEIG, EigValues, Perm )
515
         IF( MINVAL( Perm ) < 1 .OR. MAXVAL( Perm ) > NEIG ) THEN
516
           CALL Fatal(Caller,'Reordering of EigenValues failed')
517
         END IF
518

519
!
520
!        Extract the values to Elmer structures:
521
!        -----------------------------------------
522
         CALL Info( Caller, ' ', Level=4 )
523
         CALL Info( Caller, 'Eigen system solution complete: ', Level=4 )
524
         CALL Info( Caller, ' ', Level=4 )
525
         WRITE( Message,'(A,ES12.3)') 'Convergence criterion is: ', TOL
526
         CALL Info( Caller, Message, Level=7 )
527
         CALL Info( Caller,'Number of converged Ritz values is: '//I2S(IPARAM(5)),Level=4)
528
         CALL Info( Caller,'Number of update iterations taken: '//I2S(IPARAM(3)),Level=4)
529
         CALL Info( Caller,'Computed '//I2S(NEIG)//' Eigen Values',Level=4)
530

531
         ! Restore matrix values, if modified when using shift:
532
         ! ---------------------------------------------------
533
         IF ( SigmaR /= 0.0d0 ) THEN
534
           Matrix % Values = Matrix % Values + SigmaR * Matrix % MassValues
535
         END IF
536

537
         ! extract vectors:
538
         ! ----------------
539
         k = 1
540
         DO i=1,NEIG
541
           p = Perm(i)
542
           WRITE( Message,'(I0,A,2ES15.6)') i,': ',EigValues(i)
543
           CALL Info( Caller, Message, Level=4 )
544

545
           k = 1
546
           DO j=1,p-1
547
              IF ( D(j,2) == 0 ) THEN
548
                 k = k + 1
549
              ELSE
550
                 k = k + 2
551
              END IF
552
           END DO
553

554
           CALL Info(Caller,'Copying Eigenvectors to solution',Level=12)
555

556
           DO j=1,N
557
              IF ( D(p,2) /= 0.0d0 ) THEN
558
                 EigVectors(i,j) = CMPLX( V(j,k),V(j,k+1),KIND=dp )
559
              ELSE
560
                 EigVectors(i,j) = CMPLX( V(j,k),0.0d0,KIND=dp )
561
              END IF
562
           END DO
563

564
           ! Normalization moved to ScaleEigenVectors
565

566
         END DO
567
                    
568
         IF ( ListGetLogical( Params, 'Eigen System Compute Residuals', stat ) ) THEN
569
           CALL Info(Caller,'Computing eigen system residuals',Level=8)
570
           CALL CheckResiduals( Matrix, Neig, EigValues, EigVectors )
571
         END IF
572
         CALL Info( Caller, '--------------------------------',Level=4 )
573
      END IF
574

575
      DEALLOCATE( WORKL, D, WORKEV, V, CHOOSE, Perm )
576

577
#else
578
      CALL Fatal( Caller, 'Arpack Eigen System Solver not available!' )
579
#endif
580
!
581
!------------------------------------------------------------------------------
582
    END SUBROUTINE ArpackEigenSolve
583
!------------------------------------------------------------------------------
584

585

586
!------------------------------------------------------------------------------
587
!> Scale both real and complex valued eigenvectors. 
588
!------------------------------------------------------------------------------
589
    SUBROUTINE ScaleEigenVectors( Matrix, EigVectors, NoEigen, NormalizeToUnity)
590

591
      USE Multigrid
592

593
      IMPLICIT NONE
594

595
      TYPE(Matrix_t), TARGET :: Matrix
596
      COMPLEX(KIND=dp) :: EigVectors(:,:)
597
      INTEGER :: n, NoEigen
598
      LOGICAL :: NormalizeToUnity
599

600
      INTEGER :: i,j,k,l, mk, mj
601
      REAL(KIND=dp) :: r
602
      COMPLEX(KIND=dp) :: s, s1, mx
603

604
      CHARACTER(*), PARAMETER :: Caller = 'ScaleEigenVectors'
605

606
      
607
      CALL Info(Caller,'Scaling eigen vectors',Level=10)
608

609
      ! Real case: Normalize eigenvector  (x) so that x^T(M x) = 1
610
      ! Complex case: Normalize eigenvectors (x) so that x^H(M x) = 1
611
      ! (probably already done, but no harm in redoing!)
612
      ! Optionally normalize such that the maximum amplitude is set one.
613
      ! -----------------------------------------------------------------------------
614
      
615
      n = Matrix % NumberOfRows
616
      IF ( Matrix % Complex ) n = n / 2
617

618
      DO i = 1, NoEigen
619

620
        IF(  Matrix % COMPLEX ) THEN
621
          s = 0.0_dp
622
          IF( NormalizeToUnity ) THEN
623
            DO j=1,n
624
              s1 = EigVectors(i,j) * CONJG(EigVectors(i,j))
625
              IF( ABS( s1 ) > ABS( s ) ) s = s1
626
            END DO
627
          ELSE IF ( Matrix % Lumped ) THEN
628
            DO j=1,n
629
              s = s + ABS( EigVectors(i,j) )**2 * Matrix % MassValues( Matrix % Diag(2*j-1) )
630
            END DO
631
          ELSE
632
            DO j=1,Matrix % NumberOfRows,2
633
              DO l=Matrix % Rows(j),Matrix % Rows(j+1)-1,2
634
                mx = CMPLX( Matrix % MassValues(l), -Matrix % MassValues(l+1), KIND=DP )
635
                mj  = (j-1)/2 + 1
636
                mk  = (Matrix % Cols(l)-1)/2 + 1
637
                s = s + mx * CONJG( EigVectors(i,mj) ) * EigVectors(i,mk)
638
              END DO
639
            END DO
640
          END IF
641
          
642
          s = CMPLX( ParallelReduction( REAL(s) ), ParallelReduction( AIMAG(s) ), KIND=dp )
643
                  
644
          IF ( ABS(s) > 0 ) THEN
645
            s = SQRT(s) 
646
            WRITE(Message,'(A,2ES12.3)') 'Normalizing Eigenvector with: ',REAL(s),AIMAG(s)
647
            CALL Info(Caller,Message,Level=12)
648
            EigVectors(i,1:n) = EigVectors(i,1:n) / s
649
          ELSE
650
            CALL Warn(Caller,'Eigenmode has zero amplitude!')
651
          END IF
652
        ELSE          
653
          r = 0.0_dp
654
          IF( NormalizeToUnity ) THEN
655
            DO j=1,n
656
              r = MAX( r, ABS(EigVectors(i,j))**2 )
657
            END DO
658
          ELSE IF ( Matrix % Lumped ) THEN
659
            DO j=1,n
660
              r= r + ABS(EigVectors(i,j))**2 * &
661
                  Matrix % MassValues(Matrix % Diag(j))
662
            END DO
663
          ELSE
664
            r = 0
665
            DO j=1,n
666
              DO l=Matrix % Rows(j), Matrix % Rows(j+1)-1
667
                r = r +  CONJG(EigVectors(i,j)) * Matrix % MassValues(l) * EigVectors(i,Matrix % Cols(l))
668
              END DO
669
            END DO
670
          END IF
671
          
672
          r = ParallelReduction(r) 
673

674
          IF( ABS(r - 1) < EPSILON( r ) ) THEN
675
            CALL Info(Caller,'Eigenmode already normalized!',Level=12)              
676
          ELSE IF ( ABS(r) > 0 ) THEN            
677
            r = SQRT( r ) 
678
            WRITE(Message,'(A,ES12.3)') 'Normalizing Eigenvector with: ',r
679
            CALL Info(Caller,Message,Level=12)
680
            EigVectors(i,:) = EigVectors(i,:) / r
681
          ELSE
682
            CALL Warn(Caller,'Eigenmode has zero amplitude!')
683
          END IF
684
        END IF
685
          
686
      END DO
687

688
    END SUBROUTINE ScaleEigenVectors
689
!------------------------------------------------------------------------------
690

691

692
!------------------------------------------------------------------------------
693
!> Expand complex valued eigenvector to a real that is actually the same vector ;-)
694
!------------------------------------------------------------------------------
695
    SUBROUTINE ExpandEigenVectors( Matrix, EigVectors, NoEigen, dofs )
696

697
      USE Types
698

699
      IMPLICIT NONE
700

701
      TYPE(Matrix_t), TARGET :: Matrix
702
      COMPLEX(KIND=dp) :: EigVectors(:,:)
703
      INTEGER :: NoEigen
704
      INTEGER :: dofs
705

706
      COMPLEX(KIND=dp), ALLOCATABLE :: EigTmp(:)
707
      INTEGER :: n,i,j,k,cdofs
708
      COMPLEX(KIND=dp) :: s
709

710
      CHARACTER(*), PARAMETER :: Caller = 'ExpandEigenVectors'
711

712
      IF ( .NOT. Matrix % COMPLEX ) RETURN
713
    
714
      CALL Info(Caller,'Expanding eigen vectors',Level=10)
715
      
716
      n = Matrix % NumberOfRows / dofs
717
      cdofs = dofs / 2
718

719
      ALLOCATE(EigTmp(SIZE(EigVectors(1,:))))
720
      
721
      DO i = 1, NoEigen
722

723
        EigTmp = EigVectors(i,:)
724

725
        DO j = 1, n
726
          DO k = 1, cdofs          
727
            s = EigTmp(cdofs*(j-1)+k)
728
            EigVectors(i,dofs*(j-1)+k) = REAL(s)
729
            EigVectors(i,dofs*(j-1)+k+cdofs) = AIMAG(s)
730
          END DO
731
        END DO
732
      END DO
733

734
    END SUBROUTINE ExpandEigenVectors
735
!------------------------------------------------------------------------------
736

737

738
    
739
!------------------------------------------------------------------------------
740
    SUBROUTINE CheckResiduals( Matrix, n, Eigs, EigVectors )
741
!------------------------------------------------------------------------------
742
      USE CRSMatrix
743
      TYPE(Matrix_t), POINTER :: Matrix
744
      INTEGER :: i,n,sz
745
      COMPLEX(KIND=dp) :: Eigs(:), EigVectors(:,:)
746

747
      REAL(KIND=dp), ALLOCATABLE :: x(:), y(:)
748

749
      sz = Matrix % NumberOfRows
750
      ALLOCATE( x(sz), y(sz) )
751
      DO i=1,n
752
        Matrix % Values = Matrix % Values - REAL(Eigs(i)) * Matrix % MassValues
753
        x = REAL( EigVectors(i,1:sz) )
754
        CALL CRS_MatrixVectorMultiply( Matrix, x, y )
755
        Matrix % Values = Matrix % Values + REAL(Eigs(i)) * Matrix % MassValues
756

757
        WRITE( Message, * ) 'L^2 Norm of the residual: ', i, SQRT(SUM(y**2))
758
        CALL Info( 'CheckResiduals', Message, Level = 3 )
759
      END DO
760
      DEALLOCATE( x,y )
761
!------------------------------------------------------------------------------
762
END SUBROUTINE CheckResiduals
763
!------------------------------------------------------------------------------
764

765

766
!------------------------------------------------------------------------------
767
!> Solution of Eigen value problems using ARPACK library, stabilized version. 
768
!------------------------------------------------------------------------------
769
     SUBROUTINE ArpackStabEigenSolve( Solver, &
770
          Matrix, N, NEIG, EigValues, EigVectors )
771
!------------------------------------------------------------------------------
772
      USE Multigrid
773

774
      IMPLICIT NONE
775

776
      TYPE(Matrix_t), POINTER :: Matrix
777
      TYPE(Solver_t), TARGET :: Solver
778
      INTEGER :: N, NEIG
779
      INTEGER, ALLOCATABLE :: DPERM(:)
780
      COMPLEX(KIND=dp) :: EigValues(:), EigVectors(:,:)
781

782
#ifdef USE_ARPACK
783
!
784
!     %--------------%
785
!     | Local Arrays |
786
!     %--------------%
787
!
788
      REAL(KIND=dp), TARGET, ALLOCATABLE :: WORKD(:), RESID(:)
789
      REAL(KIND=dp), POINTER CONTIG :: x(:), b(:)
790
      INTEGER :: IPARAM(11), IPNTR(14)
791
      INTEGER, ALLOCATABLE :: Perm(:)
792
      LOGICAL, ALLOCATABLE :: Choose(:)
793
      REAL(KIND=dp), ALLOCATABLE :: WORKL(:), D(:,:), V(:,:)
794
!
795
!     %---------------%
796
!     | Local Scalars |
797
!     %---------------%
798
!
799
      CHARACTER ::     BMAT*1, Which*2, DirectMethod*200
800
      INTEGER   ::     IDO, NCV, lWORKL, kinfo, i, j, k, l, p, IERR, iter, &
801
                       NCONV, maxitr, ishfts, mode, istat
802
      LOGICAL   ::     First, Stat, Direct = .FALSE., FactorizeSetting, &
803
                       Iterative = .FALSE., NewSystem, &
804
                       rvec = .TRUE.
805
      LOGICAL :: Factorize, FreeFactorize,FoundFactorize,FoundFreeFactorize
806
      REAL(KIND=dp) :: SigmaR, SigmaI, TOL
807

808
      COMPLEX(KIND=dp) :: s
809
!
810
      REAL(KIND=dp), POINTER CONTIG :: SaveValues(:)
811

812
      CHARACTER(*), PARAMETER :: Caller = 'StabEigenSolve'
813

814
      
815
!     %-----------------------%
816
!     | Executable Statements |
817
!     %-----------------------%
818
!
819
!     %----------------------------------------------------%
820
!     | The number N is the dimension of the matrix. A     |
821
!     | generalized eigenvalue problem is solved (BMAT =   |
822
!     | 'G'.) NEV is the number of eigenvalues to be       |
823
!     | approximated.  The user can modify NEV, NCV, WHICH |
824
!     | to solve problems of different sizes, and to get   |
825
!     | different parts of the spectrum.  However, The     |
826
!     | following conditions must be satisfied:            |
827
!     |                     N <= MAXN,                     | 
828
!     |                   NEV <= MAXNEV,                   |
829
!     |               NEV + 1 <= NCV <= MAXNCV             | 
830
!     %----------------------------------------------------%
831
!
832
      IF ( Matrix % Lumped ) THEN
833
         CALL Error(Caller,'Lumped matrices are not allowed in stability analysis.' )
834
      END IF
835

836
      NCV = 3 * NEIG + 1
837

838
      lWORKL = NCV*(NCV+8)
839

840
      ALLOCATE(workd(3*n), resid(n), dperm(n), STAT=istat)
841
      IF ( istat /= 0 ) CALL Fatal(Caller, 'Memory allocation error.' )
842

843
      ALLOCATE( WORKL(lWORKL), D(NCV,2), V(N,NCV), CHOOSE(NCV), STAT=istat )
844
      IF ( istat /= 0 ) CALL Fatal(Caller, 'Memory allocation error.' )
845

846
      TOL = ListGetConstReal( Solver % Values, 'Eigen System Convergence Tolerance', stat )
847
      IF ( .NOT. stat ) THEN
848
         TOL = 100 * ListGetConstReal( Solver % Values, 'Linear System Convergence Tolerance' )
849
      END IF
850
!
851
!     %---------------------------------------------------%
852
!     | This program uses exact shifts with respect to    |
853
!     | the current Hessenberg matrix (IPARAM(1) = 1).    |
854
!     | IPARAM(3) specifies the maximum number of Arnoldi |
855
!     | iterations allowed.  Mode 2 of DSAUPD is used     |
856
!     | (IPARAM(7) = 2).  All these options may be        |
857
!     | changed by the user. For details, see the         |
858
!     | documentation in DSAUPD.                          |
859
!     %---------------------------------------------------%
860
!
861
      IDO   = 0
862
      kinfo = 0
863
      ishfts = 1
864
      BMAT  = 'G'
865
      Mode = 2
866

867
      SELECT CASE( ListGetString( Solver % Values,'Eigen System Select', stat ) )
868
      CASE( 'smallest magnitude' )
869
         Which = 'LM'
870
      CASE( 'largest magnitude')
871
         Which = 'SM'
872
      CASE( 'smallest real part')
873
         Which = 'LR'
874
      CASE( 'largest real part')
875
         Which = 'SR'
876
      CASE( 'smallest imag part' )
877
         Which = 'LI'
878
      CASE( 'largest imag part' )
879
         Which = 'SI'
880
      CASE( 'smallest algebraic' )
881
         Which = 'LA'
882
      CASE( 'largest algebraic' )
883
         Which = 'SA'
884
      CASE DEFAULT
885
         Which = 'LM'
886
      END SELECT
887

888
      Maxitr = ListGetInteger( Solver % Values, 'Eigen System Max Iterations', stat )
889
      IF ( .NOT. stat ) Maxitr = 300
890

891
      IPARAM = 0
892
      IPARAM(1) = ishfts
893
      IPARAM(3) = maxitr 
894
      IPARAM(7) = mode
895

896
      SigmaR = 0.0d0
897
      SigmaI = 0.0d0
898
      V = 0.0d0
899
      D = 0.0d0
900

901
      Factorize = ListGetLogical( Solver % Values, &
902
            'Linear System Refactorize', FoundFactorize )
903
      CALL ListAddLogical( Solver % Values, 'Linear System Refactorize',.TRUE. )
904

905
      FreeFactorize = ListGetLogical( Solver % Values, &
906
                'Linear System Refactorize', FoundFreeFactorize )
907
      CALL ListAddLogical( Solver % Values,  &
908
                     'Linear System Free Factorization',.FALSE. )
909

910
      Direct = ListGetString( Solver % Values, &
911
           'Linear System Solver', stat ) == 'direct'
912
      IF ( Direct ) THEN
913
         DirectMethod = ListGetString( Solver % Values, &
914
           'Linear System Direct Method', stat )
915

916
         SELECT CASE( DirectMethod )
917
         CASE('umfpack', 'big umfpack','mumps', 'superlu', 'pardiso', 'cholmod' )
918
         CASE DEFAULT
919
            Stat = CRS_ILUT(Matrix, 0.0d0)
920
         END SELECT
921
      END IF
922

923
      Iterative = ListGetString( Solver % Values, &
924
               'Linear System Solver', stat ) == 'iterative'
925

926
      stat = ListGetLogical( Solver % Values, 'No Precondition Recompute', stat  )
927

928
      IF ( Iterative .AND. Stat ) THEN
929
         CALL ListAddLogical( Solver % Values, 'No Precondition Recompute', .FALSE. )
930
      END IF
931
!
932
!     %-------------------------------------------%
933
!     | M A I N   L O O P (Reverse communication) |
934
!     %-------------------------------------------%
935
!
936
      iter = 1
937
      NewSystem = .TRUE.
938

939
      DO WHILE( ido /= 99 )
940

941
         CALL DSAUPD( ido, BMAT, n, Which, NEIG, TOL, &
942
              RESID, NCV, V, n, IPARAM, IPNTR, WORKD, WORKL, lWORKL, kinfo )
943

944
         IF( ido==-1 .OR. ido==1 ) THEN
945
            IF(InfoActive(20)) THEN
946
              CALL Info(Caller,'Arpack reverse communication calls: '//I2S(iter))
947
            ELSE
948
              CALL Info(Caller, '.', .TRUE.)
949
            END IF
950
            iter = iter + 1
951
         END IF
952

953
         SELECT CASE( ido )
954
         CASE( -1, 1 )
955
            SaveValues => Matrix % Values
956
            Matrix % Values => Matrix % MassValues
957
            CALL CRS_MatrixVectorMultiply( Matrix, WORKD(IPNTR(1)), WORKD(IPNTR(2)) )
958
            Matrix % Values => SaveValues
959
            
960
            DO i=0,n-1
961
               WORKD( IPNTR(1)+i ) = WORKD( IPNTR(2)+i )
962
            END DO
963
            
964
            IF ( Direct ) THEN
965
              CALL DirectSolver( Matrix,WORKD(IPNTR(2)),WORKD(IPNTR(1)), Solver )
966
            ELSE               
967
               x => workd(ipntr(2):ipntr(2)+n-1)
968
               b => workd(ipntr(1):ipntr(1)+n-1)
969
               
970
               IF ( Solver % MultiGridSolver ) THEN
971
                  CALL MultiGridSolve( Matrix, x, b, Solver % Variable % DOFs,  &
972
                       Solver, Solver % MultiGridLevel, NewSystem )
973
               ELSE
974
                  CALL IterSolver( Matrix, x, b, Solver )
975
               END IF
976
               
977
            END IF
978

979
         CASE( 2 )            
980
            CALL CRS_MatrixVectorMultiply( Matrix, WORKD(IPNTR(1)), WORKD(IPNTR(2)) )
981

982
         END SELECT
983

984
         IF ( NewSystem .AND. ido /= 2 ) THEN
985
            IF ( Iterative ) THEN
986
               CALL ListAddLogical( Solver % Values,  'No Precondition Recompute', .TRUE. )
987
            ELSE
988
               CALL ListAddLogical( Solver % Values, 'Linear System Refactorize', .FALSE. )
989
            END IF
990
            NewSystem = .FALSE.
991
         END IF
992

993
!-----------------------------------------------------------------------------------------      
994
      END DO  ! ido == 99
995

996
      IF ( FoundFactorize ) THEN
997
        CALL ListAddLogical( Solver % Values, 'Linear System Refactorize', Factorize )
998
      ELSE
999
        CALL ListRemove( Solver % Values, 'Linear System Refactorize' )
1000
      END IF
1001

1002
      IF ( .NOT. FoundFreeFactorize ) THEN
1003
        CALL ListRemove( Solver % Values, 'Linear System Free Factorization' )
1004
      ELSE
1005
        CALL ListAddLogical( Solver % Values, 'Linear System Free Factorization', FreeFactorize )
1006
      END IF
1007
!
1008
!     %-----------------------------------------%
1009
!     | Either we have convergence, or there is |
1010
!     | an error.                               |
1011
!     %-----------------------------------------%
1012
!
1013
      IF ( kinfo /= 0 ) THEN
1014
         WRITE( Message, * ) 'Error with DSAUPD, info = ',kinfo
1015
         CALL Fatal( 'StabEigenSolve', Message )
1016
!
1017
      ELSE 
1018
         D = 0.0d0
1019
         rvec = .TRUE.
1020

1021
         CALL DSEUPD ( rvec, 'A', Choose, D, V, N, SigmaR, &
1022
              BMAT, n, Which, NEIG, TOL, RESID, NCV, V, N, &
1023
              IPARAM, IPNTR, WORKD, WORKL, lWORKL, IERR )
1024
            
1025
!        %----------------------------------------------%
1026
!        | Eigenvalues are returned in the First column |
1027
!        | of the two dimensional array D and the       |
1028
!        | corresponding eigenvectors are returned in   |
1029
!        | the First NEV columns of the two dimensional |
1030
!        | array V if requested.  Otherwise, an         |
1031
!        | orthogonal basis for the invariant subspace  |
1032
!        | corresponding to the eigenvalues in D is     |
1033
!        | returned in V.                               |
1034
!        %----------------------------------------------%
1035
!
1036
         IF (IERR /= 0) THEN
1037
            WRITE( Message, * ) ' Error with DSEUPD, info = ', IERR
1038
            CALL Fatal( Caller, Message )
1039
         END IF
1040
!
1041
!        %------------------------------------------%
1042
!        | Print additional convergence information |
1043
!        %------------------------------------------%
1044
!
1045
         IF ( kinfo == 1 ) THEN
1046
            CALL Fatal( Caller, 'Maximum number of iterations reached.' )
1047
         ELSE IF ( kinfo == 3 ) THEN
1048
            CALL Fatal( Caller, 'No shifts could be applied during implicit Arnoldi update, try increasing NCV.' )
1049
         END IF      
1050
!
1051
!        Sort the eigenvalues to ascending order:
1052
!        ----------------------------------------
1053
         ALLOCATE( Perm(NEIG) )
1054
         Perm = [ (i, i=1,NEIG) ]
1055
         DO i=1,NEIG
1056
            EigValues(i) = CMPLX( 1.0d0 / D(i,1), D(i,2),KIND=dp )
1057
         END DO
1058

1059
         CALL SortC( NEIG, EigValues, Perm )
1060
!
1061
!        Extract the values to ELMER structures:
1062
!        -----------------------------------------
1063
         CALL Info( Caller, ' ', Level=4 )
1064
         CALL Info( Caller, 'EIGEN SYSTEM SOLUTION COMPLETE: ', Level=4 )
1065
         CALL Info( Caller, ' ', Level=4 )
1066
         WRITE( Message,'(A,ES12.3)') 'Convergence criterion is: ', TOL
1067
         CALL Info( Caller, Message, Level=7 )
1068
         CALL Info( Caller,'Number of converged Ritz values is: '//I2S(IPARAM(5)),Level=4)
1069
         CALL Info( Caller, ' ', Level=7 )
1070
         CALL Info( Caller, 'Computed Eigen Values: ', Level=4 )
1071
         CALL Info( Caller, '--------------------------------', Level=7 )
1072
         k = 1
1073

1074
         DO i=1,NEIG
1075
           p = Perm(i)
1076
           WRITE( Message, * ) i,EigValues(i)
1077
           CALL Info( Caller, Message, Level=4 )
1078

1079
           k = 1
1080
           DO j=1,p-1
1081
              IF ( D(j,2) == 0 ) THEN
1082
                 k = k + 1
1083
              ELSE
1084
                 k = k + 2
1085
              END IF
1086
           END DO
1087

1088
           DO j=1,N
1089
              IF ( D(p,2) /= 0.0d0 ) THEN
1090
                 EigVectors(i,j) = CMPLX( V(j,k),V(j,k+1),KIND=dp )
1091
              ELSE
1092
                 EigVectors(i,j) = CMPLX( V(j,k),0.0d0,KIND=dp )
1093
              END IF
1094
           END DO
1095
           
1096
           ! Normalizatin moved to ScaleEigenVectors
1097
        END DO
1098

1099
        CALL Info( Caller, '--------------------------------',Level=4 )
1100

1101
      END IF
1102

1103
      DO i = 1,Neig
1104
         IF( REAL(EigValues(i)) < 0.0d0 ) THEN
1105
            EigVectors(i,:) = EigVectors(i,:) * CMPLX(0.0d0,1.0d0,KIND=dp)
1106
         END IF
1107
      END DO
1108

1109
      DEALLOCATE( WORKL, D, V, CHOOSE, Perm )
1110
#else
1111
      CALL Fatal( Caller, 'Arpack Eigen System Solver not available!' )
1112
#endif
1113
!
1114
!------------------------------------------------------------------------------
1115
    END SUBROUTINE ArpackStabEigenSolve
1116
!------------------------------------------------------------------------------
1117

1118

1119
!------------------------------------------------------------------------------
1120
     SUBROUTINE ArpackEigenSolveComplex( Solver,Matrix,N,NEIG, &
1121
                      EigValues, EigVectors )
1122
!------------------------------------------------------------------------------
1123
!> Solution of Eigen value problems using ARPACK library, complex valued version. 
1124
!------------------------------------------------------------------------------
1125
      USE Multigrid
1126

1127
      IMPLICIT NONE
1128

1129
      TYPE(Matrix_t), POINTER :: Matrix, A
1130
      TYPE(Solver_t), TARGET :: Solver
1131
      INTEGER :: N, NEIG
1132
      COMPLEX(KIND=dp) :: EigValues(:), EigVectors(:,:)
1133

1134
#ifdef USE_ARPACK
1135
!
1136
!     %--------------%
1137
!     | Local Arrays |
1138
!     %--------------%
1139
!
1140
      COMPLEX(KIND=dp), ALLOCATABLE :: WORKD(:), RESID(:)
1141
      INTEGER :: IPARAM(11), IPNTR(14)
1142
      INTEGER, ALLOCATABLE :: Perm(:), dPerm(:)
1143
      LOGICAL, ALLOCATABLE :: Choose(:)
1144
      COMPLEX(KIND=dp), ALLOCATABLE :: WORKL(:), D(:), WORKEV(:), V(:,:)
1145

1146
!
1147
!     %---------------%
1148
!     | Local Scalars |
1149
!     %---------------%
1150
!
1151
      CHARACTER ::     BMAT*1, Which*2
1152
      INTEGER   ::     IDO, NCV, lWORKL, kinfo, i, j, k, l, p, IERR, iter, &
1153
                       NCONV, maxitr, ishfts, mode, istat, dofs
1154
      LOGICAL   ::     First, Stat, Direct = .FALSE., FoundFactorize,&
1155
                       Iterative = .FALSE., NewSystem, Factorize, FreeFactorize, FoundFreeFactorize
1156

1157
      CHARACTER(:), ALLOCATABLE :: DirectMethod, Method
1158
      COMPLEX(KIND=dp) :: Sigma = 0.0d0, s
1159
      REAL(KIND=dp), TARGET, ALLOCATABLE :: x(:), b(:), rwork(:)
1160
      REAL(KIND=dp) :: SigmaR, SigmaI, TOL
1161
!
1162
      REAL(KIND=dp), POINTER CONTIG :: SaveValues(:), SaveRhs(:)
1163

1164
      TYPE(ValueList_t), POINTER :: Params
1165

1166
      CHARACTER(*), PARAMETER :: Caller = 'EigenSolveComplex'
1167
    
1168
!
1169
!     %-----------------------%
1170
!     | Executable Statements |
1171
!     %-----------------------%
1172
!
1173
!     %----------------------------------------------------%
1174
!     | The number N is the dimension of the matrix. A     |
1175
!     | generalized eigenvalue problem is solved (BMAT =   |
1176
!     | 'G'.) NEV is the number of eigenvalues to be       |
1177
!     | approximated.  The user can modify NEV, NCV, WHICH |
1178
!     | to solve problems of different sizes, and to get   |
1179
!     | different parts of the spectrum.  However, The     |
1180
!     | following conditions must be satisfied:            |
1181
!     |                     N <= MAXN,                     | 
1182
!     |                   NEV <= MAXNEV,                   |
1183
!     |               NEV + 1 <= NCV <= MAXNCV             | 
1184
!     %----------------------------------------------------%
1185
!
1186

1187
      CALL Info(Caller,'Starting eigen system solution with complex matrices!',Level=6)
1188

1189
      Params => Solver % Values
1190

1191
      NCV = ListGetInteger( Params, 'Eigen System Lanczos Vectors', stat )
1192
      IF ( .NOT. stat ) NCV = 3*NEIG + 1
1193

1194
      IF ( NCV <=  NEIG ) THEN
1195
        CALL Fatal( 'EigenSolve', & 
1196
            'Number of Lanczos vectors must exceed the number of eigenvalues.' )
1197
      END IF
1198

1199
      ALLOCATE( workd(3*n), resid(n), dperm(n), x(2*n), b(2*n), rwork(n), STAT=istat)
1200

1201
      IF ( istat /= 0 ) THEN
1202
         CALL Fatal(Caller, 'Memory allocation error.' )
1203
      END IF
1204

1205
      ALLOCATE( WORKL(3*NCV**2 + 6*NCV), D(NCV), &
1206
         WORKEV(3*NCV), V(n,NCV+1), CHOOSE(NCV), STAT=istat )
1207

1208
      IF ( istat /= 0 ) THEN
1209
         CALL Fatal(Caller, 'Memory allocation error.' )
1210
      END IF
1211
!
1212
!     %--------------------------------------------------%
1213
!     | The work array WORKL is used in DSAUPD as        |
1214
!     | workspace.  Its dimension LWORKL is set as       |
1215
!     | illustrated below.  The parameter TOL determines |
1216
!     | the stopping criterion.  If TOL<=0, machine      |
1217
!     | precision is used.  The variable IDO is used for |
1218
!     | reverse communication and is initially set to 0. |
1219
!     | Setting INFO=0 indicates that a random vector is |
1220
!     | generated in DSAUPD to start the Arnoldi         |
1221
!     | iteration.                                       |
1222
!     %--------------------------------------------------%
1223
!
1224
!
1225
      TOL = ListGetConstReal( Solver % Values, 'Eigen System Convergence Tolerance', stat )
1226
      IF ( .NOT. stat ) THEN
1227
         TOL = 100 * ListGetConstReal( Solver % Values, 'Linear System Convergence Tolerance' )
1228
      END IF
1229

1230
      lWORKL = 3*NCV**2 + 6*NCV 
1231
      IDO   = 0
1232
      kinfo = 0
1233
!
1234
!     %---------------------------------------------------%
1235
!     | This program uses exact shifts with respect to    |
1236
!     | the current Hessenberg matrix (IPARAM(1) = 1).    |
1237
!     | IPARAM(3) specifies the maximum number of Arnoldi |
1238
!     | iterations allowed.  Mode 2 of DSAUPD is used     |
1239
!     | (IPARAM(7) = 2).  All these options may be        |
1240
!     | changed by the user. For details, see the         |
1241
!     | documentation in DSAUPD.                          |
1242
!     %---------------------------------------------------%
1243
!
1244
      ishfts = 1
1245
      BMAT  = 'G'
1246
      IF ( Matrix % Lumped ) THEN
1247
         Mode  =  2
1248
         SELECT CASE(ListGetString( Solver % Values, 'Eigen System Select',stat) )
1249
         CASE( 'smallest magnitude' )
1250
              Which = 'SM'
1251
         CASE( 'largest magnitude')
1252
              Which = 'LM'
1253
         CASE( 'smallest real part')
1254
              Which = 'SR'
1255
         CASE( 'largest real part')
1256
              Which = 'LR'
1257
         CASE( 'smallest imag part' )
1258
              Which = 'SI'
1259
         CASE( 'largest imag part' )
1260
              Which = 'LI'
1261
         CASE DEFAULT
1262
              Which = 'SM'
1263
         END SELECT
1264
      ELSE
1265
         Mode  = 3
1266
         SELECT CASE(ListGetString( Solver % Values, 'Eigen System Select',stat) )
1267
         CASE( 'smallest magnitude' )
1268
              Which = 'LM'
1269
         CASE( 'largest magnitude')
1270
              Which = 'SM'
1271
         CASE( 'smallest real part')
1272
              Which = 'LR'
1273
         CASE( 'largest real part')
1274
              Which = 'SR'
1275
         CASE( 'smallest imag part' )
1276
              Which = 'LI'
1277
         CASE( 'largest imag part' )
1278
              Which = 'SI'
1279
         CASE DEFAULT
1280
              Which = 'LM'
1281
         END SELECT
1282
      END IF
1283
!
1284
      Maxitr = ListGetInteger( Solver % Values, 'Eigen System Max Iterations', stat )
1285
      IF ( .NOT. stat ) Maxitr = 300
1286

1287
      IPARAM = 0
1288
      IPARAM(1) = ishfts
1289
      IPARAM(3) = maxitr 
1290
      IPARAM(7) = mode
1291

1292
      SigmaR = 0
1293
      SigmaI = 0
1294
      V = 0
1295

1296
!     Compute LU-factors for (A-\sigma M) (if consistent mass matrix):
1297
!     ----------------------------------------------------------------
1298
      Factorize = ListGetLogical( Params, &
1299
            'Linear System Refactorize', FoundFactorize )
1300
      CALL ListAddLogical( Params, 'Linear System Refactorize',.TRUE. )
1301

1302
      FreeFactorize = ListGetLogical( Params, &
1303
                'Linear System Refactorize', FoundFreeFactorize )
1304

1305
      CALL ListAddLogical( Params,  &
1306
                     'Linear System Free Factorization',.FALSE. )
1307

1308
      IF ( .NOT. Matrix % Lumped ) THEN
1309
        SigmaR = ListGetConstReal( Params,'Eigen System Shift', stat )
1310
        SigmaI = ListGetConstReal( Params,'Eigen System Shift Im', stat )
1311
        Sigma = CMPLX(SigmaR,SigmaI, KIND=dp)
1312

1313
        IF ( Sigma /= 0._dp ) THEN
1314
          Matrix % Values = Matrix % Values - Sigma * Matrix % MassValues
1315
        END IF
1316

1317
        Method = ListGetString( Params,'Linear System Solver', stat )         
1318
        IF ( Method == 'direct' ) THEN
1319
          DirectMethod = ListGetString( Params, &
1320
              'Linear System Direct Method', stat )
1321
          
1322
          SELECT CASE( DirectMethod )
1323
          CASE('umfpack', 'big umfpack', 'mumps', 'superlu', 'pardiso', 'cholmod')
1324
          CASE DEFAULT
1325
            Stat = CRS_ComplexILUT(Matrix, 0._dp)
1326
          END SELECT
1327
        END IF
1328
      END IF
1329

1330

1331
!
1332
!     %-------------------------------------------%
1333
!     | M A I N   L O O P (Reverse communication) |
1334
!     %-------------------------------------------%
1335
!
1336
      iter = 1
1337
      NewSystem = .TRUE.
1338

1339
      Iterative = ListGetString( Solver % Values, &
1340
        'Linear System Solver', stat ) == 'iterative'
1341

1342
      stat = ListGetLogical( Solver % Values,  'No Precondition Recompute', stat  )
1343

1344
      IF ( Iterative .AND. Stat ) THEN
1345
         CALL ListAddLogical( Solver % Values, 'No Precondition Recompute', .FALSE. )
1346
      END IF
1347

1348
      A => Matrix
1349

1350
      DO WHILE( ido /= 99 )
1351
!        %---------------------------------------------%
1352
!        | Repeatedly call the routine DSAUPD and take | 
1353
!        | actions indicated by parameter IDO until    |
1354
!        | either convergence is indicated or maxitr   |
1355
!        | has been exceeded.                          |
1356
!        %---------------------------------------------%
1357
!
1358
         CALL ZNAUPD ( ido, BMAT, n, Which, NEIG, TOL, &
1359
           RESID, NCV, v, n, IPARAM, IPNTR, WORKD, WORKL, lWORKL, RWORK, kinfo )
1360

1361
         IF (ido == -1 .OR. ido == 1) THEN
1362
            IF(InfoActive(20)) THEN
1363
              CALL Info(Caller,'Arpack reverse communication calls: '//I2S(iter))
1364
            ELSE
1365
              CALL Info(Caller, '.', .TRUE.)
1366
            END IF
1367
            Iter = Iter + 1
1368
!---------------------------------------------------------------------
1369
!           Perform  y <--- OP*x = inv[M]*A*x   (lumped mass)
1370
!                    ido =-1 inv(A-sigmaR*M)*M*x 
1371
!                    ido = 1 inv(A-sigmaR*M)*z
1372
!---------------------------------------------------------------------
1373

1374
            ! Some strategies (such as 'block') may depend on that these are set properly 
1375
            ! to reflect the linear problem under study.            
1376
            SaveRhs => A % rhs
1377
            A % rhs => b
1378
            Dofs = Solver % Variable % Dofs
1379

1380
            IF ( ido == -1 ) THEN
1381
              SaveValues => A % Values
1382
              A % Values => A % MassValues
1383
              CALL CRS_ComplexMatrixVectorMultiply( A, WORKD(IPNTR(1)), WORKD(IPNTR(2)) )
1384
              A % Values => SaveValues
1385
              
1386
              DO i=0,n-1
1387
                b(2*i+1) = REAL(  WORKD( IPNTR(2)+i ) )
1388
                b(2*i+2) = AIMAG( WORKD( IPNTR(2)+i ) )
1389
              END DO
1390
            ELSE
1391
              DO i=0,n-1
1392
                b(2*i+1) = REAL(  WORKD( IPNTR(3)+i ) )
1393
                b(2*i+2) = AIMAG( WORKD( IPNTR(3)+i ) )
1394
              END DO
1395
            END IF
1396

1397
            DO i=0,n-1
1398
               x(2*i+1) = REAL(  WORKD( IPNTR(2)+i ) )
1399
               x(2*i+2) = AIMAG( WORKD( IPNTR(2)+i ) )
1400
            END DO
1401

1402
            SELECT CASE( Method ) 
1403
            CASE('multigrid')
1404
              CALL MultiGridSolve( A, x, b, &
1405
                  DOFs, Solver, Solver % MultiGridLevel, NewSystem )
1406
            CASE('iterative')
1407
              CALL IterSolver( A, x, b, Solver )
1408
            CASE('block')
1409
              CALL BlockSolveExt( A, x, b, Solver )
1410
            CASE ('direct')
1411
              CALL DirectSolver( A, x, b, Solver )
1412
            CASE DEFAULT
1413
              CALL Fatal(Caller,'Unknown linear system method: '//TRIM(Method))
1414
            END SELECT
1415

1416
            A % rhs => SaveRhs
1417

1418
            DO i=0,n-1
1419
              WORKD( IPNTR(2)+i ) = CMPLX( x(2*i+1), x(2*i+2),KIND=dp )
1420
            END DO
1421
!
1422
         ELSE IF (ido == 2) THEN
1423
!
1424
!           %-----------------------------------------%
1425
!           |         Perform  y <--- M*x.            |
1426
!           | Need the matrix vector multiplication   |
1427
!           | routine here that takes WORKD(IPNTR(1)) |
1428
!           | as the input and returns the result to  |
1429
!           | WORKD(IPNTR(2)).                        |
1430
!           %-----------------------------------------%
1431
!
1432
            SaveValues => A % Values
1433
            A % Values => A % MassValues
1434
            CALL CRS_ComplexMatrixVectorMultiply(A, WORKD(IPNTR(1)), WORKD(IPNTR(2)) )
1435
            A % Values => SaveValues
1436
         END IF 
1437

1438
         IF ( NewSystem .AND. ido /= 2 ) THEN
1439
            IF ( Iterative ) THEN
1440
               CALL ListAddLogical( Solver % Values,  'No Precondition Recompute', .TRUE. )
1441
            ELSE
1442
               CALL ListAddLogical( Solver % Values, 'Linear System Refactorize', .FALSE. )
1443
            END IF
1444
            NewSystem = .FALSE.
1445
         END IF
1446
       END DO
1447

1448
       CALL ListAddLogical( Solver % Values, 'Linear System Refactorize', .TRUE. )
1449
       CALL ListAddLogical( Solver % Values, &
1450
                           'Linear System Free Factorization', .TRUE. )
1451
!
1452
!     %-----------------------------------------%
1453
!     | Either we have convergence, or there is |
1454
!     | an error.                               |
1455
!     %-----------------------------------------%
1456
!
1457
      IF ( kinfo /= 0 ) THEN
1458
!
1459
!        %--------------------------%
1460
!        | Error message, check the |
1461
!        | documentation in DNAUPD  |
1462
!        %--------------------------%
1463
!
1464
         WRITE( Message, * ) 'Error with DNAUPD, info = ',kinfo
1465
         CALL Fatal( Caller, Message )
1466
!
1467
      ELSE 
1468
!
1469
!        %-------------------------------------------%
1470
!        | No fatal errors occurred.                 |
1471
!        | Post-Process using DSEUPD.                |
1472
!        |                                           |
1473
!        | Computed eigenvalues may be extracted.    |  
1474
!        |                                           |
1475
!        | Eigenvectors may also be computed now if  |
1476
!        | desired.  (indicated by rvec = .true.)    | 
1477
!        %-------------------------------------------%
1478
!           
1479
         D = 0.0d0
1480
         CALL ZNEUPD ( .TRUE., 'A', Choose, D, V, N, Sigma, WORKEV, BMAT, N, Which, NEIG, &
1481
           TOL, RESID, NCV, V, N, IPARAM, IPNTR, WORKD, WORKL, lWORKL, RWORK, IERR )
1482

1483
!        %----------------------------------------------%
1484
!        | Eigenvalues are returned in the First column |
1485
!        | of the two dimensional array D and the       |
1486
!        | corresponding eigenvectors are returned in   |
1487
!        | the First NEV columns of the two dimensional |
1488
!        | array V if requested.  Otherwise, an         |
1489
!        | orthogonal basis for the invariant subspace  |
1490
!        | corresponding to the eigenvalues in D is     |
1491
!        | returned in V.                               |
1492
!        %----------------------------------------------%
1493
!
1494
         IF (IERR /= 0) THEN 
1495
!
1496
!           %------------------------------------%
1497
!           | Error condition:                   |
1498
!           | Check the documentation of DNEUPD. |
1499
!           %------------------------------------%
1500
! 
1501
            WRITE( Message, * ) ' Error with DNEUPD, info = ', IERR
1502
            CALL Fatal( Caller, Message )
1503
         END IF
1504
!
1505
!        %------------------------------------------%
1506
!        | Print additional convergence information |
1507
!        %------------------------------------------%
1508
!
1509
         IF ( kinfo == 1 ) THEN
1510
            CALL Fatal( Caller, 'Maximum number of iterations reached.' )
1511
         ELSE IF ( kinfo == 3 ) THEN
1512
            CALL Fatal( Caller, 'No shifts could be applied during implicit Arnoldi update, try increasing NCV.' )
1513
         END IF      
1514
!
1515
!        Sort the eigenvalues to ascending order:
1516
!        ----------------------------------------
1517
         ALLOCATE( Perm(NEIG) )
1518
         Perm = [ (i, i=1,NEIG) ]
1519
         DO i=1,NEIG
1520
            EigValues(i) = D(i)
1521
         END DO
1522
         CALL SortC( NEIG, EigValues, Perm )
1523

1524
!
1525
!        Extract the values to ELMER structures:
1526
!        -----------------------------------------
1527
         CALL Info( Caller, ' ', Level=4 )
1528
         CALL Info( Caller, 'EIGEN SYSTEM SOLUTION COMPLETE: ', Level=4 )
1529
         CALL Info( Caller, ' ', Level=4 )
1530
         WRITE( Message,'(A,ES12.3)') 'Convergence criterion is: ', TOL
1531
         CALL Info( Caller, Message, Level=7 )
1532
         CALL Info( Caller,'Number of converged Ritz values is: '//I2S(IPARAM(5)),Level=4)
1533
         CALL Info( Caller, ' ', Level=7 )
1534
         CALL Info( Caller, 'Computed Eigen Values: ', Level=4 )
1535
         CALL Info( Caller, '--------------------------------', Level=7 )
1536

1537
         ! Restore matrix values, if modified when using shift:
1538
         ! ---------------------------------------------------
1539
         IF ( Sigma /= 0._dp ) THEN
1540
           Matrix % Values = Matrix % Values + Sigma * Matrix % MassValues
1541
         END IF
1542

1543
         EigVectors = -1.0_dp
1544

1545
         !PRINT *,'NEIG:',NEIG,n,SIZE(EigVectors,1),SIZE(EigVectors,2),&
1546
         !    SIZE(V,1),SIZE(V,2)
1547
         
1548
         k = 1
1549
         DO i=1,NEIG
1550
            p = Perm(i)
1551
            WRITE( Message, * ) i,EigValues(i)
1552
            CALL Info( Caller, Message, Level=4 )
1553

1554
            DO j=1,N
1555
              EigVectors(i,j) = V(j,p)
1556
            END DO
1557
         END DO
1558

1559
         IF ( ListGetLogical( Params, 'Eigen System Compute Residuals', stat ) ) THEN
1560
           CALL Info(Caller,'Computing eigen system residuals',Level=8)
1561
           CALL CheckResidualsComplex( Matrix, Neig, EigValues, EigVectors )
1562
         END IF
1563
         CALL Info( Caller, '--------------------------------',Level=4 )
1564
!
1565
      END IF
1566

1567
      DEALLOCATE( WORKL, D, WORKEV, V, CHOOSE, Perm )
1568
      
1569
      CALL Info(Caller,'Finished eigen system solution!',Level=8)
1570
      
1571
#else
1572
      CALL Fatal( Caller, 'Arpack Eigen System Solver not available!' )
1573
#endif
1574
!
1575
!------------------------------------------------------------------------------
1576
     END SUBROUTINE ArpackEigenSolveComplex
1577
!------------------------------------------------------------------------------
1578

1579

1580
!------------------------------------------------------------------------------
1581
    SUBROUTINE CheckResidualsComplex( Matrix, n, Eigs, EigVectors )
1582
!------------------------------------------------------------------------------
1583
      USE CRSMatrix
1584
      TYPE(Matrix_t), POINTER :: Matrix
1585
      INTEGER :: i,j,k,n,sz
1586
      COMPLEX(KIND=dp) :: Eigs(:), EigVectors(:,:)
1587

1588
      REAL(KIND=dp), ALLOCATABLE, TARGET :: vals(:)
1589
      REAL(KIND=dp), POINTER CONTIG :: svals(:)
1590
      COMPLEX(KIND=dp) :: c,m
1591
      COMPLEX(KIND=dp), ALLOCATABLE :: x(:), y(:)
1592

1593
      sz = Matrix % NumberOfRows/2
1594
      ALLOCATE( x(sz), y(sz), vals(size(matrix % values)) ); vals=0
1595
      DO i=1,n
1596
        DO j=1,sz
1597
          DO k=Matrix % Rows(2*j-1), Matrix % Rows(2*j)-1,2
1598
            c = CMPLX(Matrix % Values(k), -Matrix % Values(k+1),KIND=dp)
1599
            m = CMPLX(Matrix % MassValues(k), -Matrix% MassValues(k+1),KIND=dp)
1600
            c = c - eigs(i) * m
1601
            vals(k) = REAL(c)
1602
            vals(k+1) = -AIMAG(c)
1603
          END DO
1604
        END DO
1605

1606
        x = EigVectors(i,:)
1607
        svals => Matrix % Values
1608
        Matrix % Values => vals
1609
        CALL CRS_ComplexMatrixVectorMultiply( Matrix, x, y )
1610
        Matrix % Values => svals
1611

1612
        WRITE( Message, * ) 'L^2 Norm of the residual: ', i, SQRT(SUM(ABS(y)**2))
1613
        CALL Info( 'CheckResidualsComplex', Message, Level = 3 )
1614
      END DO
1615
      DEALLOCATE( x,y )
1616
!------------------------------------------------------------------------------
1617
END SUBROUTINE CheckResidualsComplex
1618
!------------------------------------------------------------------------------
1619

1620

1621
!------------------------------------------------------------------------------
1622
     SUBROUTINE ArpackDampedEigenSolve( Solver, KMatrix, N, NEIG, EigValues, &
1623
          EigVectors )
1624
!------------------------------------------------------------------------------
1625
!> Solution of Eigen value problems using ARPACK library, damped version. 
1626
!------------------------------------------------------------------------------
1627
      USE Multigrid
1628
      USE ElementUtils
1629

1630
      IMPLICIT NONE
1631

1632
      TYPE(Matrix_t), POINTER :: KMatrix
1633
      TYPE(Solver_t), TARGET :: Solver
1634
      INTEGER :: N, NEIG
1635
      COMPLEX(KIND=dp) :: EigValues(:), EigVectors(:,:)
1636

1637
#ifdef USE_ARPACK
1638

1639
!     %--------------%
1640
!     | Local Arrays |
1641
!     %--------------%
1642

1643
      TYPE(Matrix_t), POINTER :: MMatrix, BMatrix
1644
      REAL(KIND=dp), POINTER CONTIG :: x(:), b(:)
1645
      INTEGER :: IPARAM(11), IPNTR(14)
1646
      REAL(KIND=dp), ALLOCATABLE, TARGET :: WORKD(:), RESID(:)
1647
      INTEGER, ALLOCATABLE :: Perm(:), kMap(:), dPerm(:)
1648
      LOGICAL, ALLOCATABLE :: Choose(:)
1649
      CHARACTER(:), ALLOCATABLE :: str
1650
      COMPLEX(KIND=dp) :: s, EigTemp(NEIG)
1651
      REAL(KIND=dp), ALLOCATABLE :: WORKL(:), D(:,:), WORKEV(:), V(:,:)
1652

1653
!     %---------------%
1654
!     | Local Scalars |
1655
!     %---------------%
1656

1657
      CHARACTER ::     BMAT*1, Which*2
1658
      INTEGER   ::     IDO, NCV, lWORKL, kinfo, i, j, k, l, p, IERR, iter, &
1659
                       NCONV, maxitr, ishfts, mode, istat, DampedMaxIter, ILU
1660
      LOGICAL   ::     First, Stat, NewSystem, UseI = .FALSE.
1661
      REAL(KIND=dp) :: SigmaR, SigmaI, TOL, DampedTOL, IScale
1662

1663
      CHARACTER(*), PARAMETER :: Caller = 'DampedEigenSolve'
1664

1665
      
1666
!     %-------------------------------------%
1667
!     | So far only iterative solver        |
1668
!     | and non-lumped matrixes are allowed |
1669
!     %-------------------------------------%
1670

1671
      IF ( KMatrix % Lumped ) THEN
1672
         CALL Error( Caller, 'Lumped matrixes are not allowed' )
1673
      END IF
1674

1675
      IF (  ListGetString( Solver % Values, 'Linear System Solver', Stat ) &
1676
           == 'direct' ) THEN
1677
         CALL Error( Caller, 'Direct solver is not allowed' )
1678
      END IF
1679

1680
      IF ( Solver % MultiGridSolver ) THEN
1681
         CALL Error( Caller, 'MultiGrid solver is not allowed' )
1682
      END IF
1683

1684
      Stat = ListGetLogical( Solver % Values, &
1685
           'No Precondition Recompute', Stat  )
1686

1687
      IF ( Stat ) THEN
1688
         CALL ListAddLogical( Solver % Values, &
1689
              'No Precondition Recompute', .FALSE. )
1690
      END IF
1691

1692

1693
!     %-----------------------%
1694
!     | Executable Statements |
1695
!     %-----------------------%
1696

1697
!     %----------------------------------------------------%
1698
!     | The number N is the dimension of the matrix. A     |
1699
!     | generalized eigenvalue problem is solved (BMAT =   |
1700
!     | 'G'.) NEV is the number of eigenvalues to be       |
1701
!     | approximated.  The user can modify NEV, NCV, WHICH |
1702
!     | to solve problems of different sizes, and to get   |
1703
!     | different parts of the spectrum.  However, The     |
1704
!     | following conditions must be satisfied:            |
1705
!     |                     N <= MAXN,                     | 
1706
!     |                   NEV <= MAXNEV,                   |
1707
!     |               NEV + 1 <= NCV <= MAXNCV             | 
1708
!     %----------------------------------------------------%
1709

1710
      NCV = 3 * NEIG + 1
1711

1712
      ALLOCATE( workd(3*n), resid(n), dperm(n), stat=istat )
1713
      IF ( istat /= 0 ) CALL Fatal( Caller, 'Memory allocation error.' )
1714

1715
      ALLOCATE( WORKL(3*NCV**2 + 6*NCV), D(NCV,3), &
1716
         WORKEV(3*NCV), V(n,NCV), CHOOSE(NCV), STAT=istat )
1717
      IF ( istat /= 0 ) CALL Fatal( Caller, 'Memory allocation error.' )
1718

1719
      CHOOSE = .FALSE.
1720
      Workl = 0.0d0
1721
      workev = 0.0d0
1722
      v = 0.0d0
1723
      d = 0.0d0
1724

1725

1726
!     %--------------------------------------------------%
1727
!     | The work array WORKL is used in DSAUPD as        |
1728
!     | workspace.  Its dimension LWORKL is set as       |
1729
!     | illustrated below.  The parameter TOL determines |
1730
!     | the stopping criterion.  If TOL<=0, machine      |
1731
!     | precision is used.  The variable IDO is used for |
1732
!     | reverse communication and is initially set to 0. |
1733
!     | Setting INFO=0 indicates that a random vector is |
1734
!     | generated in DSAUPD to start the Arnoldi         |
1735
!     | iteration.                                       |
1736
!     %--------------------------------------------------%
1737

1738
      TOL = ListGetConstReal( Solver % Values, &
1739
           'Eigen System Convergence Tolerance', Stat )
1740

1741
      IF ( .NOT. Stat ) THEN
1742
         TOL = 100 * ListGetConstReal( Solver % Values, &
1743
              'Linear System Convergence Tolerance' )
1744
      END IF
1745

1746
      DampedMaxIter = ListGetInteger( Solver % Values, &
1747
           'Linear System Max Iterations', Stat, 1 )
1748

1749
      IF ( .NOT. Stat ) DampedMaxIter = 100
1750

1751
      DampedTOL = ListGetConstReal( Solver % Values, &
1752
           'Linear System Convergence Tolerance', Stat )
1753

1754
      IF ( .NOT. Stat ) DampedTOL = TOL / 100
1755

1756
      UseI = ListGetLogical( Solver % Values, &
1757
                     'Eigen System Use Identity', Stat )
1758

1759
      IF ( .NOT. Stat ) UseI = .TRUE.
1760
!      IF ( .NOT. Stat ) UseI = .FALSE.   Changed by Antti 2004-02-18
1761

1762
      IDO   = 0
1763
      kinfo = 0
1764
      lWORKL = 3*NCV**2 + 6*NCV 
1765
!     %---------------------------------------------------%
1766
!     | This program uses exact shifts with respect to    |
1767
!     | the current Hessenberg matrix (IPARAM(1) = 1).    |
1768
!     | IPARAM(3) specifies the maximum number of Arnoldi |
1769
!     | iterations allowed.  Mode 2 of DSAUPD is used     |
1770
!     | (IPARAM(7) = 2).  All these options may be        |
1771
!     | changed by the user. For details, see the         |
1772
!     | documentation in DSAUPD.                          |
1773
!     %---------------------------------------------------%
1774
      
1775
      ishfts = 1
1776
      BMAT  = 'G'
1777
      Mode  = 3
1778
      
1779
      SELECT CASE( ListGetString(Solver % Values, 'Eigen System Select',Stat) )
1780
         CASE( 'smallest magnitude' )
1781
         Which = 'LM'
1782
         
1783
         CASE( 'largest magnitude')
1784
         Which = 'SM'
1785
         
1786
         CASE( 'smallest real part')
1787
         Which = 'LR'
1788
         
1789
         CASE( 'largest real part')
1790
         Which = 'SR'
1791
         
1792
         CASE( 'smallest imag part' )
1793
         Which = 'LI'
1794

1795
         CASE( 'largest imag part' )
1796
         Which = 'SI'
1797
         
1798
         CASE DEFAULT
1799
         Which = 'LM'
1800
      END SELECT
1801

1802
      Maxitr = ListGetInteger(Solver % Values,'Eigen System Max Iterations',Stat)
1803
      IF ( .NOT. Stat ) Maxitr = 300
1804

1805
      IPARAM = 0
1806
      IPARAM(1) = ishfts
1807
      IPARAM(3) = maxitr 
1808
      IPARAM(7) = mode
1809

1810
      SigmaR = 0.0d0
1811
      SigmaI = 0.0d0
1812
      V = 0.0d0
1813
!------------------------------------------------------------------------------
1814
! Create M and B matrixes
1815
!------------------------------------------------------------------------------
1816
      MMatrix => AllocateMatrix()
1817
      MMatrix % Values => KMatrix % MassValues
1818
      MMatrix % NumberOfRows = KMatrix % NumberOfRows
1819
      MMatrix % Rows => KMatrix % Rows
1820
      MMatrix % Cols => KMatrix % Cols
1821
      MMatrix % Diag => KMatrix % Diag
1822

1823
      IScale = MAXVAL( ABS( MMatrix % Values ) )
1824

1825
      BMatrix => AllocateMatrix()
1826
      BMatrix % NumberOfRows = KMatrix % NumberOfRows
1827
      BMatrix % Rows => KMatrix % Rows
1828
      BMatrix % Cols => KMatrix % Cols
1829
      BMatrix % Diag => KMatrix % Diag
1830
      BMatrix % Values => KMatrix % DampValues
1831
!------------------------------------------------------------------------------
1832
!     ILU Preconditioning
1833
!------------------------------------------------------------------------------
1834
      str = ListGetString( Solver % Values, 'Linear System Preconditioning', Stat )
1835

1836
      ILU = 0
1837
      IF ( .NOT. Stat ) THEN
1838
         CALL Warn( Caller, 'Using ILU0 preconditioning' )
1839
      ELSE
1840
         IF ( str == 'none' .OR. str == 'diagonal' .OR. &
1841
              str == 'ilut' .OR. str == 'multigrid' ) THEN
1842

1843
           ILU = 0
1844
           CALL Warn( Caller, 'Useing ILU0 preconditioning' )
1845
         ELSE IF ( SEQL(str,'ilu') ) THEN
1846
           IF(LEN(str)>=4) ILU = ICHAR(str(4:4)) - ICHAR('0')
1847
           IF ( ILU  < 0 .OR. ILU > 9 ) ILU = 0
1848
         ELSE
1849
           ILU = 0
1850
           CALL Warn( Caller,'Unknown preconditioner type, using ILU0' )
1851
         END IF
1852
      END IF
1853

1854
      KMatrix % Cholesky = ListGetLogical( Solver % Values,  &
1855
              'Linear System Symmetric ILU', Stat )
1856

1857
      Stat = CRS_IncompleteLU( KMatrix, ILU, Solver % Values )
1858
      IF ( .NOT. UseI ) Stat = CRS_IncompleteLU( MMatrix, ILU, Solver % Values )
1859

1860
!     %-------------------------------------------%
1861
!     | M A I N   L O O P (Reverse communication) |
1862
!     %-------------------------------------------%
1863

1864
      iter = 1
1865
      DO WHILE( ido /= 99 )
1866

1867
!        %---------------------------------------------%
1868
!        | Repeatedly call the routine DSAUPD and take | 
1869
!        | actions indicated by parameter IDO until    |
1870
!        | either convergence is indicated or maxitr   |
1871
!        | has been exceeded.                          |
1872
!        %---------------------------------------------%
1873
         CALL DNAUPD ( ido, BMAT, n, Which, NEIG, TOL, RESID, NCV, v, n, &
1874
                 IPARAM, IPNTR, WORKD, WORKL, lWORKL, kinfo )
1875

1876
         SELECT CASE(ido)
1877
         CASE(-1 )
1878
            !
1879
            ! ido =-1 inv(A)*M*x:
1880
            !----------------------------
1881
            x => workd( ipntr(1) : ipntr(1)+n-1 )
1882
            b => workd( ipntr(2) : ipntr(2)+n-1 )
1883

1884
            CALL EigenMGmv2( n/2, MMatrix, x, b, UseI, IScale )
1885
            DO i=1,n
1886
              x(i) = b(i)
1887
            END DO
1888
            CALL EigenBiCG( n, KMatrix, MMatrix, BMatrix, &
1889
                 b, x, DampedMaxIter, DampedTOL, UseI, IScale )
1890

1891
         CASE( 1 )
1892
            ! 
1893
            ! ido =-1 inv(A)*z:
1894
            !--------------------------
1895
            IF(InfoActive(20)) THEN
1896
              CALL Info(Caller,'Arpack reverse communication calls: '//I2S(iter))
1897
            ELSE
1898
              CALL Info(Caller, '.', .TRUE.)
1899
            END IF
1900
            iter = iter + 1
1901

1902
            x => workd( ipntr(2) : ipntr(2)+n-1 )
1903
            b => workd( ipntr(3) : ipntr(3)+n-1 )
1904

1905
            CALL EigenBiCG( n, KMatrix, MMatrix, BMatrix, &
1906
                  x, b, DampedMaxIter, DampedTOL, UseI,IScale )
1907

1908
         CASE( 2 )
1909
!           %-----------------------------------------%
1910
!           |         Perform  y <--- M*x.            |
1911
!           | Need the matrix vector multiplication   |
1912
!           | routine here that takes WORKD(IPNTR(1)) |
1913
!           | as the input and returns the result to  |
1914
!           | WORKD(IPNTR(2)).                        |
1915
!           %-----------------------------------------%
1916
            x => workd( ipntr(1): ipntr(1)+n-1 )
1917
            b => workd( ipntr(2): ipntr(2)+n-1 )
1918
            CALL EigenMGmv2( N/2, MMatrix, x, b, UseI, IScale )
1919

1920
         CASE DEFAULT
1921
         END SELECT
1922

1923
      END DO
1924

1925
!     %-----------------------------------------%
1926
!     | Either we have convergence, or there is |
1927
!     | an error.                               |
1928
!     %-----------------------------------------%
1929
!
1930
      IF ( kinfo /= 0 ) THEN
1931
!
1932
!        %--------------------------%
1933
!        | Error message, check the |
1934
!        | documentation in DNAUPD  |
1935
!        %--------------------------%
1936
!
1937
         WRITE( Message, * ) 'Error with DNAUPD, info = ',kinfo
1938
         CALL Fatal( Caller, Message )
1939
      ELSE 
1940
!        %-------------------------------------------%
1941
!        | No fatal errors occurred.                 |
1942
!        | Post-Process using DSEUPD.                |
1943
!        |                                           |
1944
!        | Computed eigenvalues may be extracted.    |  
1945
!        |                                           |
1946
!        | Eigenvectors may also be computed now if  |
1947
!        | desired.  (indicated by rvec = .true.)    | 
1948
!        %-------------------------------------------%
1949

1950
         D = 0.0d0
1951
         CALL DNEUPD ( .TRUE., 'A', Choose, D, D(1,2), V, N, SigmaR, SigmaI, &
1952
         WORKEV, BMAT, N, Which, NEIG, TOL, RESID, NCV, V, N, IPARAM, IPNTR, &
1953
                        WORKD, WORKL, lWORKL, IERR )
1954
!        %----------------------------------------------%
1955
!        | Eigenvalues are returned in the First column |
1956
!        | of the two dimensional array D and the       |
1957
!        | corresponding eigenvectors are returned in   |
1958
!        | the First NEV columns of the two dimensional |
1959
!        | array V if requested.  Otherwise, an         |
1960
!        | orthogonal basis for the invariant subspace  |
1961
!        | corresponding to the eigenvalues in D is     |
1962
!        | returned in V.                               |
1963
!        %----------------------------------------------%
1964

1965
         IF ( IERR /= 0 ) THEN 
1966
!           %------------------------------------%
1967
!           | Error condition:                   |
1968
!           | Check the documentation of DNEUPD. |
1969
!           %------------------------------------%
1970
            WRITE( Message, * ) ' Error with DNEUPD, info = ', IERR
1971
            CALL Fatal( Caller, Message )
1972
         END IF
1973

1974
!        %------------------------------------------%
1975
!        | Print additional convergence information |
1976
!        %------------------------------------------%
1977

1978
         IF ( kinfo == 1 ) THEN
1979
            CALL Fatal( Caller, 'Maximum number of iterations reached.' )
1980
         ELSE IF ( kinfo == 3 ) THEN
1981
            CALL Fatal( Caller, &
1982
                 'No shifts could be applied during implicit Arnoldi update, try increasing NCV.' )
1983
         END IF      
1984

1985
!        Sort the eigenvalues to ascending order:
1986
!        ( and keep in mind the corresponding vector )
1987
!        ---------------------------------------------
1988
         DO i = 1, NEIG
1989
            EigTemp(i) = CMPLX( D(i,1), D(i,2), KIND=dp )
1990
         END DO
1991

1992
         ALLOCATE( kMap( NEIG ) )
1993
         kMap(1) = 1
1994
         DO i = 2, NEIG
1995
            IF ( AIMAG( EigTemp(i-1) ) == 0 ) THEN
1996
               kMap(i) = kMap(i-1) + 1
1997
            ELSE IF ( EigTemp(i) == CONJG( EigTemp(i-1) ) ) THEN
1998
               kMap(i) = kMap(i-1)
1999
            ELSE
2000
               kMap(i) = kMap(i-1) + 2
2001
            END IF
2002
         END DO
2003

2004
         ALLOCATE( Perm( NEIG ) )
2005
         Perm = [ (i, i=1,NEIG) ]
2006
         CALL SortC( NEIG, EigTemp, Perm )
2007
         
2008
!        Extract the values to ELMER structures:
2009
!        -----------------------------------------
2010
         CALL Info( Caller, ' ', Level=4 )
2011
         CALL Info( Caller, 'EIGEN SYSTEM SOLUTION COMPLETE: ', Level=4 )
2012
         CALL Info( Caller, ' ', Level=4 )
2013
         WRITE( Message,'(A,ES12.3)') 'Convergence criterion is: ', TOL
2014
         CALL Info( Caller, Message, Level=7 )
2015
         CALL Info(Caller,'Number of converged Ritz values is: '//I2S(IPARAM(5)),Level=4)
2016
         CALL Info( Caller, ' ', Level=7 )
2017
         CALL Info( Caller, 'Computed Eigen Values: ', Level=4 )
2018
         CALL Info( Caller, '--------------------------------', Level=7 )
2019

2020
!------------------------------------------------------------------------------
2021
! Extracting the right eigen values and vectors.
2022
!------------------------------------------------------------------------------
2023

2024
! Take the first ones separately
2025
!------------------------------------------------------------------------------
2026
         EigValues(1) = EigTemp(1)         
2027
         
2028
         WRITE( Message, * ) 1,EigValues(1)
2029
         CALL Info( Caller, Message, Level=4 )
2030
         
2031
         p = Perm(1)
2032
         k = kMap(p)
2033
         IF( AIMAG( EigValues(1) ) == 0 ) THEN
2034
            DO j = 1, N/2
2035
               EigVectors(1,j) = CMPLX( V(j,k), 0.0d0,KIND=dp )
2036
            END DO
2037
         ELSE
2038
            DO j = 1, N/2
2039
               EigVectors(1,j) = CMPLX( V(j,k), V(j,k+1),KIND=dp )
2040
            END DO
2041
         END IF
2042

2043
! Then take the rest of requested values
2044
!------------------------------------------------------------------------------
2045
         l = 2
2046
         DO i = 2, NEIG/2
2047
            IF ( AIMAG( EigValues(i-1) ) /= 0 .AND. &
2048
                 ABS(AIMAG(EigTemp(l))) == ABS(AIMAG(EigValues(i-1)))) l=l+1
2049
            
2050
            EigValues(i) = EigTemp(l)
2051
            IF ( AIMAG( EigValues(i) ) < 0 ) THEN
2052
               EigValues(i) = CONJG( EigValues(i) )
2053
            END IF
2054
            
2055
            WRITE( Message, * ) i,EigValues(i)
2056
            CALL Info( Caller, Message, Level=4 )
2057
            
2058
            p = Perm(l)
2059
            k = kMap(p)
2060
            IF( AIMAG( EigValues(i) ) == 0 ) THEN
2061
               DO j = 1, N/2
2062
                  EigVectors(i,j) = CMPLX( V(j,k), 0.0d0,KIND=dp )
2063
               END DO
2064
            ELSE
2065
               DO j = 1, N/2
2066
                  EigVectors(i,j) = CMPLX( V(j,k), V(j,k+1),KIND=dp )
2067
               END DO
2068
            END IF
2069

2070
            l = l + 1
2071
         END DO
2072
               
2073
         DO i = 1, NEIG/2
2074
            s = 0.0d0
2075
            DO j=1,N/2
2076
               DO k = MMatrix % Rows(j), MMatrix % Rows(j+1)-1
2077
                  s = s + MMatrix % Values(k) * &
2078
                       CONJG( EigVectors(i,j) ) * EigVectors(i,MMatrix % Cols(k))
2079
               END DO
2080
            END DO
2081
            IF ( ABS(s) > 0 ) EigVectors(i,:) = EigVectors(i,:) / SQRT(s)
2082
         END DO
2083

2084
         CALL Warn(Caller,'Check that the scaling is not done twice if you call this!')
2085
         
2086
         ! Standard scaling moved to ScaleEigenVectors
2087

2088
         
2089
         CALL Info( Caller, '--------------------------------',Level=4 )
2090
      END IF
2091

2092
      DEALLOCATE( WORKL, D, WORKEV, V, CHOOSE, Perm )
2093

2094
      NULLIFY( MMatrix % Rows, MMatrix % Cols, MMatrix % Diag, MMatrix % Values )
2095
      CALL FreeMatrix( MMatrix )
2096

2097
      NULLIFY( BMatrix % Rows, BMatrix % Cols, BMatrix % Diag, BMatrix % Values )
2098
      CALL FreeMatrix( BMatrix )
2099
#else
2100
      CALL Fatal( Caller, 'Arpack Eigen System Solver not available!' )
2101
#endif
2102
!------------------------------------------------------------------------------
2103
     END SUBROUTINE ArpackDampedEigenSolve
2104
!------------------------------------------------------------------------------
2105

2106

2107

2108
!------------------------------------------------------------------------------
2109
    SUBROUTINE EigenBiCG( n, KMatrix, MMatrix, BMatrix, x, b, Rounds, TOL, UseI, IScale )
2110
!------------------------------------------------------------------------------
2111
      USE CRSMatrix
2112

2113
      TYPE(Matrix_t), POINTER :: KMatrix, MMatrix, BMatrix
2114
      INTEGER :: Rounds
2115
      REAL(KIND=dp) ::  TOL, IScale
2116
      REAL(KIND=dp) CONTIG :: x(:), b(:)
2117
      LOGICAL :: UseI
2118
!------------------------------------------------------------------------------
2119
      INTEGER :: i, n
2120
      REAL(KIND=dp) :: alpha, beta, omega, rho, oldrho, bnorm
2121
      REAL(KIND=dp), ALLOCATABLE :: r(:), Ri(:), P(:), V(:), S(:), &
2122
           T(:), T1(:), T2(:), Tmp(:)
2123
!------------------------------------------------------------------------------
2124

2125
      ALLOCATE(r(n), Ri(n), P(n), V(n), S(n), T(n), T1(n), t2(n), Tmp(n/2) )
2126

2127
      CALL EigenMGmv1( n/2, KMatrix, MMatrix, BMatrix, x, r, UseI, IScale )
2128
      r(1:n) = b(1:n) - r(1:n)
2129

2130
      Ri(1:n) = r(1:n)
2131
      P(1:n) = 0
2132
      V(1:n) = 0
2133
      omega  = 1
2134
      alpha  = 0
2135
      oldrho = 1
2136
      Tmp = 0.0d0
2137

2138
      bnorm = EigenMGdot( n,b,b )
2139

2140
      CALL Info( 'EigenBiCG', '--------------------' )
2141
      CALL Info( 'EigenBiCG', 'Begin BiCG iteration' )
2142
      CALL Info( 'EigenBiCG', '--------------------' )
2143

2144
      DO i=1,Rounds
2145
         rho = EigenMGdot( n, r, Ri )
2146
         
2147
         beta = alpha * rho / ( oldrho * omega )
2148
         P(1:n) = r(1:n) + beta * (P(1:n) - omega*V(1:n))
2149
!------------------------------------------------------------------------------
2150
         Tmp(1:n/2) = P(1:n/2)
2151

2152
         IF ( .NOT. UseI ) THEN
2153
            CALL CRS_LUSolve( n/2, MMatrix, Tmp(1:n/2) )
2154
         ELSE
2155
            Tmp(1:n/2) = Tmp(1:n/2) / IScale
2156
         END IF
2157

2158
         V(n/2+1:n) = Tmp(1:n/2)
2159

2160
         Tmp(1:n/2) = P(n/2+1:n)
2161
         CALL CRS_LUSolve( n/2, KMatrix, Tmp(1:n/2) )
2162
         V(1:n/2) = -1*Tmp(1:n/2)
2163

2164
         T1(1:n) = V(1:n)         
2165
         CALL EigenMGmv1( n/2, KMatrix, MMatrix, BMatrix, T1, V, UseI, IScale )
2166
!------------------------------------------------------------------------------
2167
         alpha = rho / EigenMGdot( n, Ri, V )         
2168
         S(1:n) = r(1:n) - alpha * V(1:n)         
2169
!------------------------------------------------------------------------------
2170
         Tmp(1:n/2) = S(1:n/2)
2171

2172
         IF ( .NOT. UseI ) THEN
2173
            CALL CRS_LUSolve( n/2, MMatrix, Tmp(1:n/2) )
2174
         ELSE
2175
            Tmp(1:n/2) = Tmp(1:n/2) / IScale
2176
         END IF
2177

2178
         T(n/2+1:n) = Tmp(1:n/2)
2179

2180
         Tmp(1:n/2) = S(n/2+1:n)
2181
         CALL CRS_LUSolve( n/2, KMatrix, Tmp(1:n/2) )
2182
         T(1:n/2) = -1*Tmp(1:n/2)
2183

2184
         T2(1:n) = T(1:n)
2185
         CALL EigenMGmv1( n/2, KMatrix, MMatrix, BMatrix, T2, T, UseI, IScale )
2186
!------------------------------------------------------------------------------
2187
         omega = EigenMGdot( n,T,S ) / EigenMGdot( n,T,T )         
2188
         oldrho = rho
2189
         r(1:n) = S(1:n) - omega*T(1:n)
2190
         x(1:n) = x(1:n) + alpha*T1(1:n) + omega*T2(1:n)
2191
!------------------------------------------------------------------------------
2192
         WRITE(*,*) i,EigenMGdot( n,r,r ) / bnorm
2193

2194
         IF ( EigenMGdot( n,r,r ) / bnorm < TOL ) THEN
2195
            CALL EigenMGmv1( n/2, KMatrix, MMatrix, BMatrix, x, r, UseI, IScale )
2196
            r(1:n) = b(1:n) - r(1:n)
2197

2198
            WRITE( Message,* ) 'Correct residual:', EigenMGdot( n,r,r ) / bnorm
2199
            CALL Info( 'EigenBiCG', Message )
2200

2201
            IF ( EigenMGdot( n,r,r ) / bnorm < TOL ) EXIT
2202
         END IF
2203
      END DO
2204

2205
      IF ( EigenMGdot( n,r,r ) / bnorm >= TOL ) THEN
2206
         CALL Fatal( 'EigenBiCG', 'Failed to converge' )
2207
      END IF
2208
!------------------------------------------------------------------------------
2209
    END SUBROUTINE EigenBiCG
2210
!------------------------------------------------------------------------------
2211

2212

2213

2214
!------------------------------------------------------------------------------
2215
    FUNCTION EigenMGdot( n, x, y ) RESULT(s)
2216
!------------------------------------------------------------------------------
2217
      USE Types
2218
      INTEGER :: n
2219
      REAL(KIND=dp) :: s, x(:), y(:)
2220
      
2221
      s = DOT_PRODUCT( x(1:n), y(1:n) )
2222
!------------------------------------------------------------------------------
2223
    END FUNCTION EigenMGdot
2224
!------------------------------------------------------------------------------
2225

2226

2227

2228
!------------------------------------------------------------------------------
2229
    SUBROUTINE EigenMGmv1( n, KMatrix, MMatrix, BMatrix, x, b, UseI, IScale )
2230
!------------------------------------------------------------------------------
2231
      USE CRSMatrix
2232

2233
      INTEGER :: n
2234
      TYPE(Matrix_t), POINTER :: KMatrix, MMatrix, BMatrix
2235
      REAL(KIND=dp) :: IScale
2236
      REAL(KIND=dp) CONTIG :: x(:), b(:)
2237
      LOGICAL :: UseI
2238

2239
      REAL(KIND=dp), ALLOCATABLE :: Tmp(:)
2240
       
2241
      ALLOCATE(Tmp(n))
2242

2243
      Tmp = 0.0d0
2244
      b = 0.0d0
2245

2246
      IF ( .NOT. UseI ) THEN
2247
         CALL CRS_MatrixVectorMultiply( MMatrix, x(n+1:2*n), Tmp(1:n) )
2248
         b(1:n) = b(1:n) + Tmp(1:n)
2249
      ELSE
2250
         b(1:n) = x(n+1:2*n) * IScale
2251
      END IF
2252

2253
      CALL CRS_MatrixVectorMultiply( KMatrix, x(1:n), Tmp(1:n) )
2254
      b(n+1:2*n) = b(n+1:2*n) - Tmp(1:n)
2255

2256
      CALL CRS_MatrixVectorMultiply( BMatrix, x(n+1:2*n), Tmp(1:n) )
2257
      b(n+1:2*n) = b(n+1:2*n) - Tmp(1:n)
2258
!------------------------------------------------------------------------------
2259
    END SUBROUTINE EigenMGmv1
2260
!------------------------------------------------------------------------------
2261

2262
!------------------------------------------------------------------------------
2263
    SUBROUTINE EigenMGmv2( n, MMatrix, x, b, UseI, IScale )
2264
!------------------------------------------------------------------------------
2265
      USE CRSMatrix
2266

2267
      INTEGER :: n
2268
      REAL(KIND=dp) CONTIG :: x(:), b(:)
2269
      REAL(KIND=dp) :: IScale
2270
      TYPE(Matrix_t), POINTER :: MMatrix
2271
      LOGICAL :: UseI
2272

2273
      IF ( .NOT. UseI ) THEN
2274
         CALL CRS_MatrixVectorMultiply( MMatrix, x(1:n), b(1:n) )
2275
      ELSE
2276
         b(1:n) = x(1:n) * IScale
2277
      END IF
2278
      CALL CRS_MatrixVectorMultiply( MMatrix, x(n+1:2*n), b(n+1:2*n) )
2279
!------------------------------------------------------------------------------
2280
    END SUBROUTINE EigenMGmv2
2281
!------------------------------------------------------------------------------
2282

2283
!------------------------------------------------------------------------------
2284
END MODULE EigenSolve
2285
!------------------------------------------------------------------------------
2286

2287
!> \} ElmerLib
2288

2289
Product

Resources

Company
