CoCalc -- CRSMatrix.F90

GitHub Repository: ElmerCSC/elmerfem
Path: blob/devel/fem/src/CRSMatrix.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
! *  Authors: Juha Ruokolainen
27
! *  Email:   [email protected]
28
! *  Web:     http://www.csc.fi/elmer
29
! *  Address: CSC - IT Center for Science Ltd.
30
! *           Keilaranta 14
31
! *           02101 Espoo, Finland 
32
! *
33
! *  Original Date: 01 Oct 1996
34
! *
35
! ****************************************************************************/
36

37
#include "huti_fdefs.h"
38
#include "../config.h"
39

40
!> \ingroup ElmerLib 
41
!> \{
42

43
!-----------------------------------------------------------------------------
44
!>  Module defining utility routines & matrix storage for sparse
45
!>  matrix in Compressed Row Storage (CRS) format.
46
!-----------------------------------------------------------------------------
47
MODULE CRSMatrix
48

49
  USE Lists
50

51
  IMPLICIT NONE
52

53
CONTAINS
54

55

56
!-----------------------------------------------------------------------------
57
!> Search CRS matrix for what?
58
!-----------------------------------------------------------------------------
59
  FUNCTION CRS_Search( N,Array,val ) RESULT ( Index )
60
!-----------------------------------------------------------------------------
61
    INTEGER :: N,val,Array(:)
62
!-----------------------------------------------------------------------------
63
    INTEGER :: Lower, Upper,Lou,Index
64
!-----------------------------------------------------------------------------
65
    Index = 0 
66
    Upper = N
67
    Lower = 1
68

69
    ! Handle the special case
70

71
    IF ( Upper == 0 ) RETURN
72

73
    DO WHILE( .TRUE. )
74
      IF ( Array(Lower) == val ) THEN
75
        Index = Lower
76
        EXIT
77
      ELSE IF ( Array(Upper) == val ) THEN
78
        Index = Upper
79
        EXIT
80
      END IF
81

82
      IF ( (Upper-Lower)>1 ) THEN
83
        Lou = ISHFT((Upper+Lower), -1)
84
        IF ( Array(Lou) < val ) THEN
85
          Lower = Lou
86
        ELSE
87
          Upper = Lou
88
        END IF
89
      ELSE
90
        EXIT
91
      END IF
92
    END DO
93
    
94
    RETURN
95

96
  END FUNCTION CRS_Search
97
!------------------------------------------------------------------------------
98

99

100

101
!------------------------------------------------------------------------------
102
!> Zero a CRS format matrix
103
!------------------------------------------------------------------------------
104
  SUBROUTINE CRS_ZeroMatrix(A)
105
!------------------------------------------------------------------------------
106
    TYPE(Matrix_t) :: A  !< Structure holding the matrix
107
!------------------------------------------------------------------------------
108
    INTEGER :: i, j
109
#ifdef _OPENMP
110
    ! First touch matrix values with similar access pattern as in sparse dgemv
111
    !$OMP PARALLEL DO SHARED(A) PRIVATE(j) DEFAULT(NONE)
112
    DO i=1,A % NumberOfRows
113
      DO j=A % Rows(i),A % Rows(i+1)-1
114
        A % Values(j) = REAL(0, dp)
115
      END DO
116
    END DO
117
    !$OMP END PARALLEL DO
118
#else
119
    A % Values = 0.0d0
120
#endif
121
  END SUBROUTINE CRS_ZeroMatrix
122
!------------------------------------------------------------------------------
123

124

125

126
!------------------------------------------------------------------------------
127
!> Zero given row from a CRS format matrix
128
!------------------------------------------------------------------------------
129
  SUBROUTINE CRS_ZeroRow( A,n )
130
!------------------------------------------------------------------------------
131
    TYPE(Matrix_t) :: A !< Structure holding the matrix
132
    INTEGER :: n        !< Row number to be zeroed
133
!------------------------------------------------------------------------------
134
 
135
    INTEGER :: i
136

137
    LOGICAL :: isMass, isDamp, EigenAnalysis, DampedAnalysis, HarmonicAnalysis, Found
138

139
    EigenAnalysis=.FALSE.; HarmonicAnalysis=.FALSE.
140
    DampedAnalysis = .FALSE.
141
    IF(ASSOCIATED(A % Solver)) THEN
142
       EigenAnalysis = A % Solver % NOFEigenValues > 0 .AND. &
143
           ListGetLogical( A % Solver % Values, 'Eigen Analysis',Found)
144

145
       DampedAnalysis =  EigenAnalysis .AND. &
146
           ListGetLogical( A % Solver % Values, 'Eigen System Damped', Found )
147

148
       HarmonicAnalysis = A % Solver % NOFEigenValues>0 .AND. &
149
          ListGetLogical( A % Solver % Values, 'Harmonic Analysis',Found)
150

151
    END IF
152

153
    isMass = .NOT.DampedAnalysis .AND. &
154
            (EigenAnalysis.OR.HarmonicAnalysis).AND.ASSOCIATED(A % MassValues)
155
    IF ( isMass ) &
156
      isMass = isMass .AND. SIZE(A % MassValues) == SIZE(A % Values)
157

158
    isDamp = (EigenAnalysis.OR.HarmonicAnalysis).AND.ASSOCIATED(A % DampValues)
159
    IF ( isDamp ) &
160
      isDamp = isDamp .AND. SIZE(A % DampValues) == SIZE(A % Values)
161

162
    DO i=A % Rows(n),A % Rows(n+1)-1
163
      A % Values(i) = 0.0_dp
164
    END DO
165

166
    IF ( isMass ) THEN
167
      DO i=A % Rows(n),A % Rows(n+1)-1
168
         A % MassValues(i) = 0.0_dp
169
      END DO
170
    END IF
171

172
    IF ( isDamp )  THEN
173
      DO i=A % Rows(n),A % Rows(n+1)-1
174
        A % DampValues(i) = 0.0_dp
175
      END DO
176
    END IF
177

178
  END SUBROUTINE CRS_ZeroRow
179
!------------------------------------------------------------------------------
180
  
181

182

183
!------------------------------------------------------------------------------
184
!>    Sort columns to ascending order for rows of a CRS format matrix-
185
!>    Optionally also sort the corresponding values of the matrix.
186
!------------------------------------------------------------------------------
187
  SUBROUTINE CRS_SortMatrix( A, ValuesToo )
188
!------------------------------------------------------------------------------
189
    TYPE(Matrix_t) :: A   !< Structure holding the matrix
190
    LOGICAL, OPTIONAL, INTENT(IN) :: ValuesToo  !< Should the values be sorted as well?
191
!------------------------------------------------------------------------------
192

193
    INTEGER :: i,j,n
194
    LOGICAL :: SortValues
195
    REAL(KIND=dp), POINTER :: Values(:)
196
    INTEGER, POINTER :: Cols(:),Rows(:),Diag(:)
197

198
    SortValues = .FALSE.
199
    IF ( PRESENT(ValuesToo) ) SortValues=ValuesToo
200

201
    Diag   => A % Diag
202
    Rows   => A % Rows
203
    Cols   => A % Cols
204
    IF ( SortValues ) Values => A % Values
205
    n = A % NumberOfRows
206

207
    IF ( .NOT. A % Ordered ) THEN
208
      IF ( SortValues ) THEN
209
        !$OMP PARALLEL DO DEFAULT(NONE) &
210
        !$OMP SHARED(Rows, Cols, Values, N) &
211
        !$OMP PRIVATE(i)
212
        DO i=1,N
213
          CALL SortF( Rows(i+1)-Rows(i),Cols(Rows(i):Rows(i+1)-1),Values(Rows(i):Rows(i+1)-1) )
214
        END DO
215
        !$OMP END PARALLEL DO
216
      ELSE
217
        !$OMP PARALLEL DO DEFAULT(NONE) &
218
        !$OMP SHARED(Rows, Cols, N) &
219
        !$OMP PRIVATE(i)
220
        DO i=1,N
221
          CALL Sort( Rows(i+1)-Rows(i),Cols(Rows(i):Rows(i+1)-1) )
222
        END DO
223
        !$OMP END PARALLEL DO
224
      END IF
225

226
      IF ( ASSOCIATED(Diag) ) THEN
227
        !$OMP PARALLEL DO DEFAULT(NONE) &
228
        !$OMP SHARED(Diag, Rows, Cols, N) &
229
        !$OMP PRIVATE(i,j)
230
        DO i=1,n
231
          DO j=Rows(i),Rows(i+1)-1
232
            IF ( Cols(j) == i ) THEN
233
              Diag(i) = j
234
              EXIT
235
            END IF
236
          END DO
237
        END DO
238
        !$OMP END PARALLEL DO
239
      END IF
240

241
      A % Ordered = .TRUE.
242
    END IF
243

244
  END SUBROUTINE CRS_SortMatrix
245
!------------------------------------------------------------------------------
246

247

248
!------------------------------------------------------------------------------
249
!>    Sort columns to ascending order for rows of a CRS format matrix-
250
!>    Optionally also sort the corresponding values of the matrix.
251
!>    This operates just on a basic matrix type.
252
!------------------------------------------------------------------------------
253
  SUBROUTINE CRS_SortBasicMatrix( A, ValuesToo )
254
!------------------------------------------------------------------------------
255
    TYPE(BasicMatrix_t), TARGET :: A  !< Structure holding the matrix
256
    LOGICAL, OPTIONAL, INTENT(IN) :: ValuesToo    !< Should the values be sorted as well?
257
!------------------------------------------------------------------------------
258

259
    INTEGER :: i,j,n
260

261
    LOGICAL :: SortValues
262
    REAL(KIND=dp), POINTER :: Values(:)
263
    INTEGER, POINTER :: Cols(:),Rows(:),Diag(:)
264

265
    SortValues = .FALSE.
266
    IF ( PRESENT(ValuesToo) ) SortValues=ValuesToo
267

268
    Diag   => A % Diag
269
    Rows   => A % Rows
270
    Cols   => A % Cols
271
    IF ( SortValues ) Values => A % Values
272
    n = A % NumberOfRows
273

274
    IF ( SortValues ) THEN
275
      DO i=1,N
276
        CALL SortF( Rows(i+1)-Rows(i),Cols(Rows(i):Rows(i+1)-1),Values(Rows(i):Rows(i+1)-1) )
277
      END DO
278
    ELSE
279
      DO i=1,N
280
          CALL Sort( Rows(i+1)-Rows(i),Cols(Rows(i):Rows(i+1)-1) )
281
      END DO
282
    END IF
283

284
    IF ( ASSOCIATED(Diag) ) THEN
285
      DO i=1,N
286
        DO j=Rows(i),Rows(i+1)-1
287
          IF ( Cols(j) == i ) THEN
288
            Diag(i) = j
289
            EXIT
290
          END IF
291
        END DO
292
      END DO
293
    END IF
294
  END SUBROUTINE CRS_SortBasicMatrix
295
!------------------------------------------------------------------------------
296

297

298

299
!------------------------------------------------------------------------------
300
!>    Fill in the column number to a CRS format matrix (values are not 
301
!>    affected in any way).
302
!------------------------------------------------------------------------------
303
  SUBROUTINE CRS_MakeMatrixIndex( A,i,j,prev )
304
!------------------------------------------------------------------------------
305
    TYPE(Matrix_t) :: A  !< Structure holding matrix
306
    INTEGER, INTENT(IN) :: i         !< row number of the matrix element
307
    INTEGER, INTENT(IN) :: j         !< column number of the matrix element
308
    INTEGER, OPTIONAL :: prev
309
!------------------------------------------------------------------------------
310

311
    INTEGER :: k,l,n
312
    INTEGER, POINTER :: Cols(:),Rows(:)
313

314
    Rows   => A % Rows
315
    Cols   => A % Cols
316

317
    n = Rows(i)
318
    l = Rows(i)
319
    IF(PRESENT(prev)) l=prev+1
320
    DO k=l,Rows(i+1)-1
321
      IF ( Cols(k) == j ) THEN
322
        RETURN
323
      ELSE IF ( Cols(k) < 1 ) THEN
324
        n = k
325
        EXIT
326
      END IF
327
    END DO
328

329
    IF ( Cols(n) >= 1 ) THEN
330
      WRITE( Message, * ) 'Trying to access non-existent column:',n,Cols(n)
331
      CALL Error( 'MakeMatrixIndex', Message )
332
      RETURN
333
    END IF
334

335
    Cols(n) = j
336
    IF(PRESENT(prev)) prev=n
337
  END SUBROUTINE CRS_MakeMatrixIndex
338
!------------------------------------------------------------------------------
339

340

341

342
!------------------------------------------------------------------------------
343
!>    Add a given value to an element of a  CRS format matrix.
344
!------------------------------------------------------------------------------
345
  SUBROUTINE CRS_AddToMatrixElement( A,i,j,val,previ )
346
!------------------------------------------------------------------------------
347
    TYPE(Matrix_t) :: A     !< Structure holding the matrix
348
    INTEGER, INTENT(IN) :: i         !< row number of the matrix element
349
    INTEGER, INTENT(IN) :: j         !< column number of the matrix element
350
    INTEGER, INTENT(INOUT), OPTIONAL :: previ     !< speed sequential access
351
    REAL(KIND=dp), INTENT(IN) :: val   !< value to be added to the matrix element
352
 !------------------------------------------------------------------------------
353
    INTEGER :: k
354
    REAL(KIND=dp), POINTER :: Values(:)
355
    INTEGER, POINTER :: Cols(:),Rows(:),Diag(:)
356
!------------------------------------------------------------------------------
357
    IF(i>A % NumberOfRows) THEN
358
      WRITE(Message,'(A,ES12.3)') 'Nonexistent row index for matrix entry:',val
359
      CALL Warn('CRS_AddToMatrixElement',Message)
360
      CALL Warn('CRS_AddToMatrixElement','Row: '//i2s(i)//' Col: '//i2s(j))
361
      CALL Warn('CRS_AddToMatrixElement','Number of Matrix rows:'//i2s(A % NumberOfRows))
362
      CALL Warn('CRS_AddToMatrixElement','Converting CRS to list')
363
      A % FORMAT = MATRIX_LIST; RETURN
364
    END IF
365

366
    Rows   => A % Rows
367
    Cols   => A % Cols
368
    Diag   => A % Diag
369
    Values => A % Values
370

371
    IF ( .NOT.ASSOCIATED(Diag) .OR. i /= j .OR. .NOT. A % Ordered ) THEN
372
      IF(PRESENT(Previ)) THEN
373
        k = previ+1
374
        DO WHILE(k<Rows(i+1)-1)
375
          IF(Cols(k)==j) EXIT
376
          k = k+1 
377
        END DO
378
        previ = k
379
        IF(Cols(k)/=j) RETURN
380
      ELSE
381
        k = CRS_Search( Rows(i+1)-Rows(i),Cols(Rows(i):Rows(i+1)-1),j )
382
        IF ( k==0 .AND. val/=0 ) THEN
383
          WRITE(Message,'(A,ES12.3)') 'Nonexistent col index for matrix entry:',val
384
          CALL Warn('CRS_AddToMatrixElement',Message)
385
          CALL Warn('CRS_AddToMatrixElement','Row: '//i2s(i)//' Col: '//i2s(j))
386
          CALL Warn('CRS_AddToMatrixElement','Number of Matrix rows:'//i2s(A % NumberOfRows))
387
          CALL Warn('CRS_AddToMatrixElement','Converting CRS to list')
388
          A % FORMAT = MATRIX_LIST
389
        END IF
390
        IF ( k==0 ) RETURN
391
        k = k + Rows(i) - 1
392
      END IF
393
    ELSE
394
      k = Diag(i)
395
    END IF
396
!$omp atomic
397
    Values(k) = Values(k) + val
398
  END SUBROUTINE CRS_AddToMatrixElement
399
!------------------------------------------------------------------------------
400

401

402
!------------------------------------------------------------------------------
403
!> Check existence of a matrix element.
404
!------------------------------------------------------------------------------
405
  FUNCTION CRS_CheckMatrixElement( A,i,j ) RESULT ( Found ) 
406
!------------------------------------------------------------------------------
407
    TYPE(Matrix_t) :: A                !< Structure holding the matrix
408
    INTEGER, INTENT(IN) :: i           !< row number of the matrix element
409
    INTEGER, INTENT(IN) :: j           !< column number of the matrix element
410
    LOGICAL :: Found                   !< Does the matrix element exist!
411
 !------------------------------------------------------------------------------
412
    INTEGER, POINTER :: Cols(:),Rows(:)
413
!------------------------------------------------------------------------------
414

415
    Found = .FALSE.
416
    IF(i > A % NumberOfRows) RETURN
417

418
    Rows => A % Rows
419
    Cols => A % Cols
420

421
    Found = ANY( Cols(Rows(i):Rows(i+1)-1) == j )
422
    
423
  END FUNCTION CRS_CheckMatrixElement
424
!------------------------------------------------------------------------------
425

426

427
  
428
!------------------------------------------------------------------------------
429
!> Check whether matrix has a symmetric topology
430
!------------------------------------------------------------------------------
431
  SUBROUTINE CRS_CheckSymmetricTopo( A )
432
!------------------------------------------------------------------------------
433
    TYPE(Matrix_t) :: A     !< Structure holding the matrix
434
!------------------------------------------------------------------------------
435
    INTEGER :: i,j,k,k2,ns
436
    INTEGER, POINTER :: Cols(:),Rows(:)
437
    LOGICAL :: Hit
438
!------------------------------------------------------------------------------
439
    Rows   => A % Rows
440
    Cols   => A % Cols
441

442
    ns = 0
443
    
444
    DO i=1,A % NumberOfRows
445
      DO k=Rows(i),Rows(i+1)-1
446
        j=Cols(k)
447
        Hit = .FALSE.
448
        DO k2=Rows(j),Rows(j+1)-1
449
          IF(Cols(k2)==i) THEN
450
            Hit = .TRUE.
451
            EXIT
452
          END IF
453
        END DO
454
        IF(.NOT. Hit) THEN
455
          ns = ns + 1
456
          !PRINT *,'Not symmetric: ',i,j
457
        END IF
458
      END DO
459
    END DO
460
    
461
    CALL Info('CSR_CheckSymmetricTopo','Number of symmetry misses:'//I2S(ns))
462
    
463
  END SUBROUTINE CRS_CheckSymmetricTopo
464
!------------------------------------------------------------------------------
465

466

467
!------------------------------------------------------------------------------
468
!> Check whether matrix has a symmetric topology
469
!------------------------------------------------------------------------------
470
  SUBROUTINE CRS_CheckComplexTopo( A )
471
!------------------------------------------------------------------------------
472
    TYPE(Matrix_t) :: A     !< Structure holding the matrix
473
!------------------------------------------------------------------------------
474
    INTEGER :: i,j,k,i2,j2,k2,nc,nr
475
    LOGICAL :: ImRow,ImCol
476
    INTEGER, POINTER :: Cols(:),Rows(:)
477
    LOGICAL :: Hit
478
!------------------------------------------------------------------------------
479
    Rows   => A % Rows
480
    Cols   => A % Cols
481

482
    nr = 0
483
    nc = 0
484
    
485
    DO i=1,A % NumberOfRows
486
      ImRow = (MODULO(i,2)==0)
487
      IF(ImRow) THEN
488
        i2=i-1
489
      ELSE
490
        i2=i+1
491
      END IF
492

493
      DO k=Rows(i),Rows(i+1)-1
494
        j=Cols(k)
495
        ImCol = (MODULO(j,2)==0) 
496
        IF(ImCol) THEN
497
          j2=j-1
498
        ELSE
499
          j2=j+1
500
        END IF
501

502
        ! We should find complementary entry on each row
503
        Hit = .FALSE.
504
        DO k2=Rows(i),Rows(i+1)-1
505
          IF(Cols(k2)==j2) THEN
506
            Hit = .TRUE.
507
            EXIT
508
          END IF
509
        END DO
510
        IF(.NOT. Hit) THEN
511
          nr = nr + 1          
512
          !PRINT *,'No complement on row: ',i,j,j2,ImRow,ImCol
513
        END IF
514

515
        ! We should find complementary entry on each column
516
        Hit = .FALSE.        
517
        DO k2=Rows(i2),Rows(i2+1)-1
518
          IF(Cols(k2)==j) THEN
519
            Hit = .TRUE.
520
            EXIT
521
          END IF
522
        END DO
523
        IF(.NOT. Hit) THEN
524
          nc = nc + 1          
525
          !PRINT *,'No complement on column: ',i,j,i2,ImRow,ImCol
526
        END IF
527
        
528
      END DO
529
    END DO
530
    
531
    CALL Info('CSR_CheckComplexTopo','Number of row misses:'//I2S(nr))
532
    CALL Info('CSR_CheckComplexTopo','Number of col misses:'//I2S(nc))
533
    
534
  END SUBROUTINE CRS_CheckComplexTopo
535
!------------------------------------------------------------------------------
536

537

538

539
  
540

541
!------------------------------------------------------------------------------
542
!>    Set a given value to an element of a  CRS format matrix.
543
!------------------------------------------------------------------------------
544
  SUBROUTINE CRS_SetMatrixElement( A,i,j,val )
545
!------------------------------------------------------------------------------
546
    TYPE(Matrix_t) :: A     !< Structure holding the matrix
547
    INTEGER, INTENT(IN) :: i         !< row number of the matrix element
548
    INTEGER, INTENT(IN) :: j         !< column number of the matrix element
549
    REAL(KIND=dp), INTENT(IN) :: val   !< new value of the matrix element
550
!------------------------------------------------------------------------------ 
551
    INTEGER :: k
552
    REAL(KIND=dp), POINTER :: Values(:)
553
    INTEGER, POINTER :: Cols(:),Rows(:),Diag(:)
554
!------------------------------------------------------------------------------
555

556
    IF(i>A % NumberOfRows) THEN
557
      A % FORMAT=MATRIX_LIST; RETURN
558
    END IF
559

560
    Rows   => A % Rows
561
    Cols   => A % Cols
562
    Diag   => A % Diag
563
    Values => A % Values
564

565
    IF ( .NOT.ASSOCIATED(Diag) .OR. i /= j .OR. .NOT. A % Ordered ) THEN
566
       k = CRS_Search( Rows(i+1)-Rows(i),Cols(Rows(i):Rows(i+1)-1),j )
567
      IF ( k==0 ) THEN
568
         A % FORMAT = MATRIX_LIST
569
         RETURN
570
       END IF
571
       k = k + Rows(i) - 1
572
    ELSE
573
       k = Diag(i)
574
    END IF
575
    Values(k) = val
576
  END SUBROUTINE CRS_SetMatrixElement
577
!------------------------------------------------------------------------------
578

579

580
!------------------------------------------------------------------------------
581
!>    Get a given matrix entry from CRS format matrix.
582
!------------------------------------------------------------------------------
583
  FUNCTION CRS_GetMatrixElement( A,i,j ) RESULT ( val )
584
!------------------------------------------------------------------------------
585
    TYPE(Matrix_t), INTENT(IN):: A     !< Structure holding the matrix
586
    INTEGER, INTENT(IN) :: i         !< row number of the matrix element
587
    INTEGER, INTENT(IN) :: j         !< column number of the matrix element
588
    REAL(KIND=dp) :: val   !< obtained value of the matrix element
589
!------------------------------------------------------------------------------ 
590
    INTEGER :: k
591
    REAL(KIND=dp), POINTER :: Values(:)
592
    INTEGER, POINTER :: Cols(:),Rows(:),Diag(:)
593
!------------------------------------------------------------------------------
594
    Rows   => A % Rows
595
    Cols   => A % Cols
596
    Diag   => A % Diag
597
    Values => A % Values
598

599
    val = REAL(0,dp)
600
    IF ( .NOT.ASSOCIATED(Diag).OR.i /= j .OR. .NOT. A % Ordered ) THEN
601
       k = CRS_Search( Rows(i+1)-Rows(i),Cols(Rows(i):Rows(i+1)-1),j )
602
      IF ( k==0 ) THEN
603
         PRINT*,'Trying to get value to nonexistent matrix element: ', i,j
604
         RETURN
605
       END IF
606
       k = k + Rows(i) - 1
607
    ELSE
608
       k = Diag(i)
609
    END IF
610
    val = Values(k)
611

612
  END FUNCTION CRS_GetMatrixElement
613
!------------------------------------------------------------------------------
614

615
!------------------------------------------------------------------------------
616
!>    Get a given matrix entry from CRS format matrix and replace it with a new value
617
!------------------------------------------------------------------------------
618
  FUNCTION CRS_ChangeMatrixElement( A,i,j, NewVal ) RESULT ( OldVal )
619
!------------------------------------------------------------------------------
620
    TYPE(Matrix_t), INTENT(IN):: A     !< Structure holding the matrix
621
    INTEGER, INTENT(IN) :: i         !< row number of the matrix element
622
    INTEGER, INTENT(IN) :: j         !< column number of the matrix element
623
    REAL(KIND=dp), INTENT(IN) :: NewVal  !< Value to be set   
624
    REAL(KIND=dp) :: OldVal !< Value to be gotten  
625
!------------------------------------------------------------------------------
626
 
627
    INTEGER :: k
628
    REAL(KIND=dp), POINTER :: Values(:)
629
    INTEGER, POINTER :: Cols(:),Rows(:),Diag(:)
630
!------------------------------------------------------------------------------
631

632
    Rows   => A % Rows
633
    Cols   => A % Cols
634
    Diag   => A % Diag
635
    Values => A % Values
636

637
    OldVal = REAL(0, dp)
638
    IF ( .NOT.ASSOCIATED(Diag).OR.i /= j .OR. .NOT. A % Ordered ) THEN
639
       k = CRS_Search( Rows(i+1)-Rows(i),Cols(Rows(i):Rows(i+1)-1),j )
640
      IF ( k==0 ) THEN
641
         PRINT*,'Trying to change value of a nonexistent matrix element: ', i,j,NewVal
642
         RETURN
643
       END IF
644
       k = k + Rows(i) - 1
645
    ELSE
646
       k = Diag(i)
647
    END IF
648
!$omp critical
649
    OldVal = Values(k)
650
    Values(k) = NewVal
651
!$omp end critical
652

653
  END FUNCTION CRS_ChangeMatrixElement
654
!------------------------------------------------------------------------------
655

656

657
!------------------------------------------------------------------------------
658
!>    Add a row together with another row of a CRS matrix, and thereafter zero it.
659
!------------------------------------------------------------------------------
660
  SUBROUTINE CRS_MoveRow( A,n1,n2,coeff,staycoeff,movecoeff )
661
!------------------------------------------------------------------------------
662
    TYPE(Matrix_t) :: A    !< Structure holding the matrix
663
    INTEGER, INTENT(IN) :: n1         !< Row number to be copied and zerod
664
    INTEGER, INTENT(IN) :: n2         !< Row number to be added
665
    REAL(KIND=dp),OPTIONAL :: coeff   !< Optional coefficient to multiply the row to be copied with
666
    REAL(KIND=dp),OPTIONAL :: staycoeff   !< Optional coefficient to multiply the staying row
667
    REAL(KIND=dp),OPTIONAL :: movecoeff   !< Optional coefficient to multiply the row copied to other direction
668
!------------------------------------------------------------------------------
669
    REAL(KIND=dp) :: val, c, d
670
    INTEGER :: i,j,i2
671
    REAL(KIND=dp), ALLOCATABLE :: Row2(:)
672
    
673
    
674
    ! memorize the row that will be written over
675
    IF( PRESENT( movecoeff ) ) THEN      
676
      i = A % Rows(n2+1)-A % Rows(n2)
677
      ALLOCATE( Row2(i) )
678
      DO i = A % Rows(n2), A % Rows(n2+1)-1
679
        i2 = i - A % Rows(n2) + 1 
680
        j = A % Cols(i)
681
        val = A % Values(i) 
682
        Row2(i2) = val
683
      END DO
684
    END IF
685

686
    
687
    IF( PRESENT(Coeff)) THEN
688
      c = coeff
689
    ELSE
690
      c = 1.0_dp
691
    END IF
692

693
    IF( PRESENT(StayCoeff)) THEN
694
      d = StayCoeff
695
    ELSE
696
      d = 0.0_dp
697
    END IF
698
    
699
    DO i=A % Rows(n1),A % Rows(n1+1)-1
700
      j = A % Cols(i)
701
      val = A % Values(i) 
702
      IF( ABS( val ) > TINY( val ) ) THEN
703
        A % Values(i) = d * val 
704
        CALL CRS_AddToMatrixElement( A,n2,j,c*val )      
705
      END IF
706
    END DO
707

708
    IF( PRESENT( movecoeff ) ) THEN      
709
      DO i = A % Rows(n2), A % Rows(n2)-1
710
        i2 = i - A % Rows(n2) + 1 
711
        j = A % Cols(i)
712
        val = Row2(i2)
713
        IF( ABS( val ) > TINY( val ) ) THEN
714
          CALL CRS_AddToMatrixElement( A,n1,j,movecoeff*val )      
715
        END IF
716
      END DO
717
    END IF
718

719
    
720
  END SUBROUTINE CRS_MoveRow
721
!------------------------------------------------------------------------------
722
  
723

724

725
!------------------------------------------------------------------------------
726
!>    Add a set of values (.i.e. element stiffness matrix) to a CRS format
727
!>    matrix. 
728
!------------------------------------------------------------------------------
729
  SUBROUTINE CRS_GlueLocalMatrix( A,N,Dofs,Indeces,LocalMatrix,GlobalValues )
730
!------------------------------------------------------------------------------
731
     TYPE(Matrix_t) :: A  !< Structure holding matrix
732
     REAL(KIND=dp), INTENT(IN) :: LocalMatrix(:,:)  !< A (N x Dofs) x ( N x Dofs) matrix holding the values to be added to the CRS format matrix
733
     INTEGER, INTENT(IN) :: N             !< Number of nodes in element
734
     INTEGER, INTENT(IN) :: Dofs          !< Number of degrees of freedom for one node
735
     INTEGER, INTENT(IN) :: Indeces(:)    !< Maps element node numbers to global (or partition) node numbers 
736
                                          !! (to matrix rows and columns, if Dofs = 1)
737
     REAL(KIND=dp), OPTIONAL, TARGET :: GlobalValues(:)
738
!------------------------------------------------------------------------------ 
739
     INTEGER :: i,j,k,l,c,Row,Col
740
     REAL(KIND=dp), POINTER :: Values(:)
741
     INTEGER, POINTER :: Cols(:),Rows(:),Diag(:)
742
!------------------------------------------------------------------------------
743

744
     Diag   => A % Diag
745
     Rows   => A % Rows
746
     Cols   => A % Cols
747
     IF(PRESENT(GlobalValues)) THEN
748
       Values => GlobalValues
749
     ELSE
750
       Values => A % Values
751
     END IF
752

753
     IF ( Dofs == 1 ) THEN
754
       DO i=1,n
755
         Row = Indeces(i)
756
         IF ( Row <=0 ) CYCLE
757
         IF ( Diag(Row) <= 0) CYCLE
758
         DO j=1,n
759
           Col = Indeces(j)
760
           IF ( Col <= 0 ) CYCLE
761
           IF (Col >= Row ) THEN
762
             DO c=Diag(Row),Rows(Row+1)-1
763
               IF ( Cols(c) == Col ) THEN
764
!$omp atomic
765
                 Values(c) = Values(c) + LocalMatrix(i,j)
766
                 EXIT
767
               END IF
768
             END DO
769
           ELSE
770
             DO c=Rows(Row),Diag(Row)-1
771
               IF ( Cols(c) == Col ) THEN
772
!$omp atomic
773
                 Values(c) = Values(c) + LocalMatrix(i,j)
774
                 EXIT
775
               END IF
776
             END DO
777
           END IF
778
         END DO
779
       END DO
780
     ELSE
781
       DO i=1,N
782
          DO k=0,Dofs-1
783
             IF ( Indeces(i) <= 0 ) CYCLE
784
             Row  = Dofs * Indeces(i) - k
785
             IF ( Diag(Row) <= 0)  CYCLE
786
             DO j=1,N
787
                DO l=0,Dofs-1
788
                   IF ( Indeces(j) <= 0 ) CYCLE
789
                   Col  = Dofs * Indeces(j) - l
790
                   IF ( Col >= Row ) THEN
791
                     DO c=Diag(Row),Rows(Row+1)-1
792
                        IF ( Cols(c) == Col ) THEN
793
!$omp atomic
794
                           Values(c) = Values(c) + LocalMatrix(Dofs*i-k,Dofs*j-l)
795
                           EXIT
796
                        END IF
797
                     END DO
798
                   ELSE
799
                     DO c=Rows(Row),Diag(Row)-1
800
                        IF ( Cols(c) == Col ) THEN
801
!$omp atomic
802
                           Values(c) = Values(c) + LocalMatrix(Dofs*i-k,Dofs*j-l)
803
                           EXIT
804
                        END IF
805
                     END DO
806
                   END IF
807
                END DO
808
             END DO
809
          END DO
810
       END DO
811
     END IF
812

813
  END SUBROUTINE CRS_GlueLocalMatrix
814
!------------------------------------------------------------------------------
815

816
  
817
  SUBROUTINE CRS_GlueLocalMatrixVec(Gmtr, N, NDOFs, Indices, Lmtr, MCAssembly, MaskedAssembly)
818
    TYPE(Matrix_t) :: Gmtr                   !< Global matrix
819
    INTEGER, INTENT(IN) :: N                 !< Number of nodes in element
820
    INTEGER, INTENT(IN) :: NDOFs             !< Number of degrees of freedom for one node
821
    INTEGER, INTENT(IN) CONTIG :: Indices(:) !< Maps element node numbers to global (or partition) node numbers
822
    REAL(KIND=dp), INTENT(IN) CONTIG :: Lmtr(:,:)   !< A (N x Dofs) x ( N x Dofs) matrix holding the values to be added
823
    LOGICAL :: MCAssembly                    !< Is the assembly multicolored or not (free of race conditions)
824
    LOGICAL :: MaskedAssembly                !< Does the assembly need masking for indices
825

826
    ! Local storage
827
    INTEGER :: Lind((N*NDOFs)*(N*NDOFs))
828
    REAL(KIND=dp) :: Lvals((N*NDOFs)*(N*NDOFs))
829
    INTEGER :: pind(N)
830

831
    INTEGER :: i,j, nzind
832
    INTEGER :: ci, ri, rli, rti, rdof, cdof, nidx, pnidx
833

834
    INTEGER, POINTER CONTIG :: gia(:), gja(:)
835
    REAL(KIND=dp), POINTER CONTIG :: gval(:)
836
!DIR$ ATTRIBUTES ALIGN:64::Lind
837
!DIR$ ATTRIBUTES ALIGN:64::Lvals
838

839
    gia  => Gmtr % Rows
840
    gja   => Gmtr % Cols
841
    gval => Gmtr % Values
842
    
843
    pnidx = -8
844
    
845
    ! Get permutation such that Indices(pind(1:N)) is sorted
846
!DIR$ INLINE
847
    CALL InsertionSort(N, Indices, pind)
848

849
    ! Check if vector masking is needed
850
    IF (MaskedAssembly) THEN
851
      ! Masking and counting needed in the assignment (slower)
852
      IF (NDOFs == 1) THEN
853
        ! Separate case for only 1 DOF per node
854

855
        ! Construct index array
856
        nzind = 0
857
        DO i=1,N
858
          IF (Indices(pind(i)) > 0) THEN
859
            ! Row index
860
            ri = Indices(pind(i))
861

862
            ! Get row pointers
863
            rli = gia(ri)
864
            rti = gia(ri+1)-1
865
            DO j=1,N
866
              ! Get global matrix index for entry (ri,Indices(j)).
867
              IF (Indices(pind(j)) > 0) THEN
868
                nzind = nzind + 1
869
                ! Get global matrix index for entry (ri,Indices(j)).
870
!DIR$ INLINE
871
                nidx = GetNextIndex(gja,Indices(pind(j)), rli, rti)
872
                Lind(nzind)=nidx
873
                Lvals(nzind)=Lmtr(pind(i),pind(j))
874
#ifdef __INTEL_COMPILER
875
                ! Issue prefetch for every new cache line of gval(nidx)
876
                IF (nidx > pnidx+8) THEN
877
                  CALL MM_PREFETCH(gval(nidx),2)
878
                  pnidx = nidx
879
                END IF
880
#endif
881
              END IF
882
            END DO
883
          END IF
884
        END DO
885
      ELSE
886
        ! More than 1 DOF per node
887
        ! Construct index array
888
        nzind = 0
889
        DO i=1,N
890
          IF (Indices(pind(i)) > 0) THEN
891
            DO rdof=1,NDOFs
892
              ! Row index
893
              ri = NDOFs*(Indices(pind(i))-1)+rdof
894

895
              ! Get row pointers
896
              rli = gia(ri)
897
              rti = gia(ri+1)-1
898
              DO j=1,N
899
                IF (Indices(pind(j)) > 0) THEN
900
                  DO cdof=1,NDOFs
901
                    ci = NDOFs*(Indices(pind(j))-1)+cdof
902
                    ! Get global matrix index for entry (ri,ci).
903
!DIR$ INLINE
904
                    nidx=GetNextIndex(gja, ci, rli, rti)
905
                    nzind = nzind + 1
906
                    Lind(nzind)=nidx
907
                    Lvals(nzind)=Lmtr(NDOFs*(pind(i)-1)+rdof,&
908
                                           NDOFs*(pind(j)-1)+cdof)
909
#ifdef __INTEL_COMPILER
910
                    ! Issue prefetch for every new cache line of gval(nidx)
911
                    IF (nidx > pnidx+8) THEN
912
                      CALL MM_PREFETCH(gval(nidx),2)
913
                      pnidx = nidx
914
                    END IF
915
#endif
916
                  END DO
917
                END IF
918
              END DO
919
            END DO
920
          END IF
921
        END DO
922
      END IF ! NDOFs==1 check
923

924
    ELSE
925
      ! No masking needed (faster)
926
      IF (NDOFs == 1) THEN
927
        ! Separate case for only 1 DOF per node
928

929
        ! Construct index array
930
        DO i=1,N
931
          ! Row index
932
          ri = Indices(pind(i))
933

934
          ! Get row pointers
935
          rli = gia(ri)
936
          rti = gia(ri+1)-1
937
          DO j=1,N
938
            ! Get global matrix index for entry (ri,Indices(j)).
939
!DIR$ INLINE
940
            nidx=GetNextIndex(gja,Indices(pind(j)), rli, rti)
941
            Lind(N*(i-1)+j)=nidx
942
            Lvals(N*(i-1)+j)=Lmtr(pind(i),pind(j))
943
#ifdef __INTEL_COMPILER
944
            ! Issue prefetch for every new cache line of gval(nidx)
945
            IF (nidx > pnidx+8) THEN
946
              CALL MM_PREFETCH(gval(nidx),2)
947
              pnidx = nidx
948
            END IF
949
#endif
950
          END DO
951
        END DO
952
        nzind = N*N
953
      ELSE
954
        ! More than 1 DOF per node
955

956
        ! Construct index array
957
        nzind = 0
958
        DO i=1,N
959
          DO rdof=1,NDOFs
960
            ! Row index
961
            ri = NDOFs*(Indices(pind(i))-1)+rdof
962

963
            ! Get row pointers
964
            rli = gia(ri)
965
            rti = gia(ri+1)-1
966
            DO j=1,N
967
              DO cdof=1,NDOFs
968
                ci = NDOFs*(Indices(pind(j))-1)+cdof
969
                ! Get global matrix index for entry (ri,ci).
970
!DIR$ INLINE
971
                nidx = GetNextIndex(gja, ci, rli, rti)
972
                nzind = nzind + 1
973
                Lind(nzind) = nidx
974
                Lvals(nzind) = Lmtr(NDOFs*(pind(i)-1)+rdof, NDOFs*(pind(j)-1)+cdof)
975
#ifdef __INTEL_COMPILER
976
                ! Issue prefetch for every new cache line of gval(nidx)
977
                IF (nidx > pnidx+8) THEN
978
                  CALL MM_PREFETCH(gval(nidx),2)
979
                  pnidx = nidx
980
                END IF
981
#endif
982
              END DO
983
            END DO
984
          END DO
985
        END DO
986
      END IF ! NDOFs==1 check
987
    END IF ! Masking check
988

989
    ! The actual contribution loop
990
    IF (MCAssembly) THEN
991
      !_ELMER_OMP_SIMD
992
      DO i=1,nzind
993
        gval(Lind(i)) = gval(Lind(i)) + Lvals(i)
994
      END DO
995
    ELSE
996
      DO i=1,nzind
997
        !$OMP ATOMIC
998
        gval(Lind(i)) = gval(Lind(i)) + Lvals(i)
999
      END DO
1000
    END IF
1001

1002
  CONTAINS
1003
    
1004
    PURE FUNCTION BinarySearch(arr, key, lind, tind) RESULT(keyloc)
1005
      IMPLICIT NONE
1006

1007
      INTEGER, INTENT(IN) CONTIG :: arr(:)
1008
      INTEGER, INTENT(IN) :: key, lind, tind
1009

1010
      INTEGER, PARAMETER :: LINSEARCHTHRESH = 8
1011
      INTEGER :: keyloc
1012
      INTEGER :: li, ti, mi
1013

1014
      li = lind
1015
      ti = tind
1016
      DO WHILE ((li+LINSEARCHTHRESH)<ti)
1017
        ! Compute midpoint
1018
        mi = li + ((ti - li) / 2)
1019

1020
        IF (arr(mi)<key) THEN
1021
          li = mi + 1
1022
        ELSE
1023
          ti = mi
1024
        END IF
1025
      END DO
1026

1027
      IF (li<ti) THEN
1028
        keyloc = 0
1029
        DO mi=li,ti
1030
          IF (arr(mi)==key) keyloc = mi
1031
        END DO
1032
      ELSE IF (li == ti .AND. arr(li)==key) THEN
1033
        keyloc = li
1034
      ELSE
1035
        keyloc = 0
1036
      END IF
1037
    END FUNCTION BinarySearch
1038
    
1039
    ! Find index matching key from arr(lind:tind). lind is set to location of 
1040
    ! arr(keyloc))=key, i.e., keyloc once the search ends
1041
    FUNCTION GetNextIndex(arr, key, lind, tind) RESULT(keyloc)
1042
      IMPLICIT NONE
1043

1044
      INTEGER, INTENT(IN) CONTIG :: arr(:)
1045
      INTEGER, INTENT(IN) :: key
1046
      INTEGER, INTENT(INOUT) :: lind
1047
      INTEGER, INTENT(IN) :: tind
1048
      INTEGER :: keyloc
1049

1050
      INTEGER :: ci
1051

1052
!DIR$ NOVECTOR
1053
      DO ci=lind,tind
1054
         IF (arr(ci)==key) EXIT
1055
      END DO
1056
      keyloc = ci
1057
      lind = keyloc
1058
    END FUNCTION GetNextIndex
1059

1060
  END SUBROUTINE CRS_GlueLocalMatrixVec
1061

1062
  SUBROUTINE InsertionSort(N, val, ind)
1063
    IMPLICIT NONE
1064
    INTEGER, INTENT(in) :: N, val(N)
1065
    INTEGER, INTENT(inout) :: ind(N)
1066
    INTEGER :: tmp, i, j
1067
    
1068
    ind(1)=1
1069
    DO i=2,N
1070
      tmp=i
1071
      ! Make room to move val(ind(i)) to its final place
1072
      ! in the sorted sequence val(ind(1)) ... val(ind(i-1))
1073
      DO j=i-1,1,-1
1074
        IF (val(ind(j))<=val(tmp)) EXIT
1075
        ind(j+1)=ind(j)
1076
      END DO
1077

1078
      ind(j+1)=tmp
1079
    END DO
1080
  END SUBROUTINE InsertionSort
1081

1082
!------------------------------------------------------------------------------
1083
!>    Add a set of values (.i.e. element stiffness matrix) to a CRS format
1084
!>    matrix. For this matrix the entries are ordered so that first for one
1085
!>    dof you got all nodes, and then for second etc. There may be an offset
1086
!>    to the entries making the subroutine suitable for coupled monolithic
1087
!>    matrix assembly.
1088
!------------------------------------------------------------------------------
1089
  SUBROUTINE CRS_GlueLocalSubMatrix( A,row0,col0,Nrow,Ncol,RowInds,ColInds,&
1090
                  RowDofs,ColDofs,LocalMatrix )
1091
!------------------------------------------------------------------------------ 
1092
     REAL(KIND=dp), INTENT(IN) :: LocalMatrix(:,:)  !< A (Nrow x RowDofs) x ( Ncol x ColDofs) matrix holding the values to be
1093
                                                    !!            added to the CRS format matrix
1094
     TYPE(Matrix_t) :: A           !< Structure holding matrix
1095
     INTEGER, INTENT(IN) :: Nrow   !< Number of nodes for the row variable
1096
     INTEGER, INTENT(IN) :: Ncol   !< Number of nodes for the column variable
1097
     INTEGER, INTENT(IN) :: RowDofs  !< Number of dofs for row variable
1098
     INTEGER, INTENT(IN) :: ColDofs  !< Number of dofs for column variable
1099
     INTEGER, INTENT(IN) :: Col0     !< Offset for column variable
1100
	 INTEGER, INTENT(IN) :: Row0     !< Offset for row variable
1101
     INTEGER, INTENT(IN) :: RowInds(:)  !< Permutation of the row dofs
1102
	 INTEGER, INTENT(IN) :: ColInds(:)  !< Permutation of the column dofs
1103
!------------------------------------------------------------------------------
1104
     INTEGER :: i,j,k,l,c,Row,Col
1105
     INTEGER, POINTER :: Cols(:),Rows(:),Diag(:)
1106
     REAL(KIND=dp), POINTER :: Values(:)
1107
!------------------------------------------------------------------------------
1108
     Diag   => A % Diag
1109
     Rows   => A % Rows
1110
     Cols   => A % Cols
1111
     Values => A % Values
1112
     
1113
     DO i=1,Nrow
1114
        DO k=0,RowDofs-1
1115
           IF ( RowInds(i) <= 0 ) CYCLE
1116
           Row  = Row0 + RowDofs * RowInds(i) - k
1117

1118
           DO j=1,Ncol
1119
              DO l=0,ColDofs-1
1120
                 IF ( ColInds(j) <= 0 ) CYCLE
1121
                 Col  = Col0 + ColDofs * ColInds(j) - l
1122

1123
! If Diag does not exist then one cannot separate the gluing into two parts
1124
! In fact this cannot be guarateed for the off-diagonal block matrices and hence 
1125
! this condition is set active.
1126
                IF( .TRUE. ) THEN
1127
                    DO c=Rows(Row),Rows(Row+1)-1
1128
                       IF ( Cols(c) == Col ) THEN
1129
!$omp atomic
1130
                           Values(c) = Values(c) + LocalMatrix(RowDofs*i-k,ColDofs*j-l)
1131
                           EXIT
1132
                        END IF
1133
                     END DO
1134
                     IF( Cols(c) /= Col ) PRINT *,'NO HIT 1',Row,Col
1135
                 ELSE IF ( Col >= Row ) THEN
1136
                    DO c=Diag(Row),Rows(Row+1)-1
1137
                       IF ( Cols(c) == Col ) THEN
1138
!$omp atomic
1139
                          Values(c) = Values(c) + LocalMatrix(RowDofs*i-k,ColDofs*j-l)
1140
                          EXIT
1141
                       END IF
1142
                    END DO
1143
                    IF( Cols(c) /= Col ) PRINT *,'NO HIT 2',Row,Col
1144
                 ELSE
1145
                    DO c=Rows(Row),Diag(Row)-1
1146
                       IF ( Cols(c) == Col ) THEN
1147
!$omp atomic
1148
                           Values(c) = Values(c) + LocalMatrix(RowDofs*i-k,ColDofs*j-l)
1149
                           EXIT
1150
                        END IF
1151
                     END DO
1152
                     IF( Cols(c) /= Col ) PRINT *,'NO HIT 3',Row,Col
1153
                  END IF
1154
               END DO
1155
            END DO
1156
         END DO
1157
      END DO
1158
    END SUBROUTINE CRS_GlueLocalSubMatrix
1159
!------------------------------------------------------------------------------
1160

1161

1162
!------------------------------------------------------------------------------
1163
!> When Dirichlet conditions are set by zeroing the row except for setting 
1164
!> the diagonal entry to one, the matrix symmetry is broken. This routine 
1165
!> maintains the symmetric structure of the matrix equation.
1166
!------------------------------------------------------------------------------
1167
  SUBROUTINE CRS_SetSymmDirichlet( A,b,n,val,s)
1168
!------------------------------------------------------------------------------
1169
    TYPE(Matrix_t) :: A       !< Structure holding matrix
1170
    INTEGER, INTENT(IN) :: n              !< Index of the dofs to be fixed   
1171
    REAL(KIND=dp) :: b(:)     !< right-hand-side of the matrix equation
1172
    REAL(KIND=dp), INTENT(IN) :: val      !< Dirichlet value to be set
1173
    REAL(KIND=dp), OPTIONAL :: s
1174
!------------------------------------------------------------------------------
1175
    INTEGER :: i,j,k,l,k1,k2
1176
    REAL(KIND=dp) :: t,ss
1177
    LOGICAL :: isMass, isDamp
1178

1179
    IF( PRESENT( s ) ) THEN
1180
      ss = s
1181
    ELSE
1182
      ss = 1.0_dp
1183
    END IF
1184
    
1185
    isMass = ASSOCIATED(A % MassValues)
1186
    IF ( isMass ) &
1187
      isMass = isMass .AND. SIZE(A % MassValues) == SIZE(A % Values)
1188

1189
    isDamp = ASSOCIATED(A % DampValues)
1190
    IF ( isDamp ) &
1191
      isDamp = isDamp .AND. SIZE(A % DampValues) == SIZE(A % Values)
1192

1193
!    IF(.NOT.A % NoDirichlet) THEN
1194
    DO l=A % Rows(n),A % Rows(n+1)-1
1195
      i = A % Cols(l)
1196

1197
      ! Cycle the diagonal entry that will be the last to set
1198
      IF ( i == n ) CYCLE
1199

1200
      ! The range of indexes to scan over assuming they are sorted
1201
      IF ( n > i ) THEN
1202
        k1 = A % Diag(i)+1
1203
        k2 = A % Rows(i+1)-1
1204
      ELSE 
1205
        k1 = A % Rows(i)
1206
        k2 = A % Diag(i)-1
1207
      END IF
1208

1209
      k = k2 - k1 + 1
1210
      IF ( k <= 30 ) THEN
1211
        DO j = k1, k2
1212
          IF ( A % Cols(j) == n ) THEN
1213
            b(i) = b(i) - A % Values(j) * val
1214
            A % Values(j) = 0.0_dp
1215
            IF ( isMass ) A % MassValues(j) = 0._dp
1216
            IF ( isDamp ) A % DampValues(j) = 0._dp
1217
            EXIT
1218
          ELSE IF ( A % Cols(j) > n ) THEN
1219
            EXIT
1220
          END IF
1221
        END DO
1222
      ELSE
1223
        j = CRS_Search( k,A % Cols(k1:k2),n )
1224
        IF ( j > 0 ) THEN
1225
          j = j + k1 - 1
1226
          b(i) = b(i) - A % Values(j) * val
1227
          A % Values(j) = 0.0_dp
1228
          IF ( isMass ) A % MassValues(j) = 0._dp
1229
          IF ( isDamp ) A % DampValues(j) = 0._dp
1230
        END IF
1231
      END IF
1232
    END DO
1233

1234
    CALL CRS_ZeroRow(A,n)
1235
    A % Values(A % Diag(n)) = ss
1236
    b(n) = ss * val
1237
       !     END IF
1238
     
1239
    !IF(ALLOCATED(A % Dvalues)) THEN
1240
      !A % DValues(n) = val
1241
    !ELSE
1242
    !  b(n) = s * val
1243
    !END IF
1244
      !IF(ALLOCATED(A % ConstrainedDOF))
1245
      !A % ConstrainedDOF(n) = .TRUE.
1246
!------------------------------------------------------------------------------
1247
  END SUBROUTINE CRS_SetSymmDirichlet
1248
!------------------------------------------------------------------------------
1249

1250
!------------------------------------------------------------------------------
1251
!> When Dirichlet conditions are set by zeroing the row except for setting 
1252
!> the diagonal entry to one, the matrix symmetry is broken. This routine 
1253
!> maintains the symmetric structure of the matrix equation.
1254
!> This routine different from the one above in that only the matrix entries
1255
!> NOT on the row of the dirichlet condition are set. 
1256
!------------------------------------------------------------------------------
1257
  SUBROUTINE CRS_ElimSymmDirichlet(A,b)
1258
!------------------------------------------------------------------------------
1259
    TYPE(Matrix_t) :: A       !< Structure holding matrix
1260
    REAL(KIND=dp) :: b(:)     !< right-hand-side of the matrix equation
1261
!------------------------------------------------------------------------------
1262
    INTEGER :: i,j,k,l,n
1263
    REAL(KIND=dp) :: t,val
1264
    LOGICAL :: isMass, isDamp
1265
    
1266
    isMass = ASSOCIATED(A % MassValues)
1267
    IF ( isMass ) &
1268
        isMass = isMass .AND. SIZE(A % MassValues) == SIZE(A % Values)
1269
    
1270
    isDamp = ASSOCIATED(A % DampValues)
1271
    IF ( isDamp ) &
1272
        isDamp = isDamp .AND. SIZE(A % DampValues) == SIZE(A % Values)
1273
    
1274
    
1275
    DO n=1,A % NumberOfRows
1276

1277
      ! There is no point eliminating entries in a row that will be nullified in the end
1278
      IF( A % ConstrainedDOF(n) ) CYCLE
1279
      
1280
      DO l=A % Rows(n),A % Rows(n+1)-1
1281
        i = A % Cols(l)
1282
        IF( A % ConstrainedDOF(i) ) THEN         
1283
          b(n) = b(n) - A % Values(l) * A % DValues(i)
1284

1285
          A % Values(l) = 0.0_dp
1286
          IF ( isMass ) A % MassValues(l) = 0._dp
1287
          IF ( isDamp ) A % DampValues(l) = 0._dp
1288
        END IF
1289
      END DO
1290
    END DO
1291

1292
    !------------------------------------------------------------------------------
1293
  END SUBROUTINE CRS_ElimSymmDirichlet
1294
!------------------------------------------------------------------------------
1295

1296
  
1297
!------------------------------------------------------------------------------
1298
!> Computes the rowsum of a given row in a CRS matrix.
1299
!------------------------------------------------------------------------------
1300
FUNCTION CRS_RowSum( A,k ) RESULT(rsum)
1301
!------------------------------------------------------------------------------
1302
   TYPE(Matrix_t), INTENT(IN) :: A       !< Structure holding matrix
1303
   INTEGER, INTENT(IN) :: k              !< Row index 
1304
   REAL(KIND=dp) :: rsum                 !< Sum of the entries on the row
1305
!------------------------------------------------------------------------------
1306
   INTEGER :: i
1307
   
1308
   rsum = 0.0D0
1309
   DO i=A % Rows(k), A % Rows(k+1)-1
1310
     rsum  = rsum + A % Values( i )
1311
   END DO
1312
!------------------------------------------------------------------------------
1313
END FUNCTION CRS_RowSum
1314
!------------------------------------------------------------------------------
1315

1316
!------------------------------------------------------------------------------
1317
!> Computes the absolute rowsum of a given row in a CRS matrix.
1318
!------------------------------------------------------------------------------
1319
FUNCTION CRS_RowSumAbs( A,k ) RESULT(rsum)
1320
!------------------------------------------------------------------------------
1321
  TYPE(Matrix_t), INTENT(IN) :: A       !< Structure holding matrix
1322
  INTEGER, INTENT(IN) :: k              !< Row index 
1323
  REAL(KIND=dp) :: rsum                 !< Sum of the entries on the row
1324
!------------------------------------------------------------------------------
1325
  INTEGER :: i
1326

1327
  rsum = 0.0D0
1328
  DO i=A % Rows(k), A % Rows(k+1)-1
1329
    rsum  = rsum + ABS( A % Values( i ) ) 
1330
  END DO
1331
!------------------------------------------------------------------------------
1332
END FUNCTION CRS_RowSumAbs
1333
!------------------------------------------------------------------------------
1334

1335

1336
!------------------------------------------------------------------------------
1337
!> Computes information on the matrix rowsums.
1338
!------------------------------------------------------------------------------
1339
SUBROUTINE CRS_RowSumInfo( A, Values ) 
1340
!------------------------------------------------------------------------------
1341
   TYPE(Matrix_t), INTENT(IN) :: A       
1342
   REAL(KIND=dp), POINTER, OPTIONAL :: Values(:)
1343
!------------------------------------------------------------------------------
1344
   REAL(KIND=dp), POINTER :: PValues(:)
1345
   INTEGER :: i,j,k              
1346
   REAL(KIND=dp) :: val,rsum,absrsum     
1347
   REAL(KIND=dp) :: minrsum,maxrsum,minabsrsum,maxabsrsum
1348
!------------------------------------------------------------------------------
1349

1350
   minrsum = HUGE(minrsum)
1351
   maxrsum = -HUGE(maxrsum)
1352
   minabsrsum = HUGE(minabsrsum)
1353
   maxabsrsum = 0.0_dp
1354
   
1355
   IF( PRESENT( Values ) ) THEN
1356
     PValues => Values
1357
   ELSE
1358
     PValues => A % Values
1359
   END IF
1360

1361
   DO i=1,A % NumberOfRows
1362
     rsum = 0.0_dp
1363
     absrsum = 0.0_dp
1364

1365
     DO j=A % Rows(i), A % Rows(i+1)-1
1366
       val = PValues( j )
1367
       rsum = rsum + val
1368
       absrsum = absrsum + ABS( val ) 
1369
     END DO
1370

1371
     minrsum = MIN( minrsum, rsum ) 
1372
     maxrsum = MAX( maxrsum, rsum ) 
1373
     minabsrsum = MIN( minabsrsum, absrsum ) 
1374
     maxabsrsum = MAX( maxabsrsum, absrsum ) 
1375
   END DO
1376

1377
   WRITE( Message,'(A,ES12.4)') 'Total sum:',SUM( PValues )
1378
   CALL Info( 'CRS_RowSumInfo', Message ) 
1379
   WRITE( Message,'(A,2ES12.4)') 'Rowsum range:',minrsum,maxrsum
1380
   CALL Info( 'CRS_RowSumInfo', Message ) 
1381
   WRITE( Message,'(A,2ES12.4)') 'Absolute rowsum range:',minabsrsum,maxabsrsum
1382
   CALL Info( 'CRS_RowSumInfo', Message ) 
1383
   
1384
!------------------------------------------------------------------------------
1385
 END SUBROUTINE CRS_RowSumInfo
1386
!------------------------------------------------------------------------------
1387

1388

1389

1390
!------------------------------------------------------------------------------
1391
!>    Create the structures required for a CRS format matrix.
1392
!------------------------------------------------------------------------------
1393
  FUNCTION CRS_CreateMatrix( N,Total,RowNonzeros,Ndeg,Reorder,AllocValues,SetRows ) RESULT(A)
1394
!------------------------------------------------------------------------------
1395
    INTEGER, INTENT(IN) :: N    !< Number of rows in the matrix
1396
    INTEGER, INTENT(IN) :: Total  !< Total number of nonzero entries in the matrix
1397
    INTEGER, INTENT(IN) :: Ndeg   !< Negrees of freedom
1398
    INTEGER, INTENT(IN), OPTIONAL :: RowNonzeros(:)  !< Number of nonzero entries in rows of the matrix
1399
    INTEGER, INTENT(IN) :: Reorder(:)      !< Permutation index for bandwidth reduction    
1400
    LOGICAL, INTENT(IN) :: AllocValues     !< Should the values arrays be allocated ?
1401
    LOGICAL, INTENT(IN), OPTIONAL :: SetRows
1402
    TYPE(Matrix_t), POINTER :: A  !>  Pointer to the created Matrix_t structure.
1403
!------------------------------------------------------------------------------
1404
    INTEGER :: i,j,k,istat
1405
    INTEGER, POINTER CONTIG :: InvPerm(:)
1406
    LOGICAL :: SetRowSizes
1407
!------------------------------------------------------------------------------
1408

1409
    CALL Info('CRS_CreateMatrix','Creating CRS Matrix of size: '//I2S(n),Level=12)
1410

1411
    SetRowSizes = .TRUE.
1412
    IF( PRESENT( SetRows ) ) SetRowSizes = SetRows
1413

1414
    A => AllocateMatrix()
1415

1416
    ALLOCATE( A % Rows(n+1),A % Diag(n),STAT=istat )
1417
    IF ( istat /= 0 ) THEN
1418
      CALL Fatal( 'CRS_CreateMatrix', 'Memory allocation error for matrix topology of size: '&
1419
          //I2S(n))
1420
    END IF
1421

1422
    k = Ndeg*Ndeg*Total
1423
    CALL Info('CRS_CreateMatrix','Creating CRS Matrix with nofs: '//I2S(k),Level=14)
1424
    ALLOCATE( A % Cols(k),STAT=istat )
1425
    IF ( istat /= 0 ) THEN
1426
      CALL Fatal( 'CRS_CreateMatrix', 'Memory allocation error for matrix cols of size: '&
1427
          //I2S(k) )
1428
    END IF
1429

1430
    IF ( AllocValues ) THEN
1431
      ALLOCATE( A % Values(k), STAT=istat )
1432
      IF ( istat /= 0 ) THEN
1433
        CALL Fatal( 'CRS_CreateMatrix', 'Memory allocation error for matrix values of size: '&
1434
            //I2S(k) )
1435
      END IF
1436
    END IF
1437

1438
    NULLIFY( A % ILUValues )
1439
    NULLIFY( A % CILUValues )
1440

1441
    A % ndeg = ndeg
1442
    A % NumberOfRows = n
1443
    A % Rows(1) = 1
1444
    A % Ordered = .FALSE.
1445

1446
    ! We don't always want to set the rows as it is more easily done elsewhere
1447
    ! but for backward compatibility the default way is maintained.
1448
    IF(.NOT. SetRowSizes ) THEN
1449
      A % Cols = 0
1450
      A % Diag = 0
1451
      RETURN
1452
    END IF
1453

1454
    
1455
    InvPerm => A % Diag ! just available memory space...
1456
    j = 0
1457
    DO i=1,SIZE(Reorder)
1458
       IF ( Reorder(i) > 0 ) THEN
1459
          j = j + 1
1460
          InvPerm(Reorder(i)) = j
1461
       END IF
1462
    END DO
1463

1464
    !$OMP PARALLEL SHARED(A, k, N, Ndeg, RowNonzeros, InvPerm) &
1465
    !$OMP PRIVATE(j) DEFAULT(NONE)
1466

1467
#ifdef _OPENMP
1468
    IF (omp_get_num_threads() > 1) THEN
1469
      ! First touch matrix values with similar access pattern as in sparse dgemv
1470
      !$OMP DO
1471
      DO i=1,N+1
1472
        A % Rows(i) = REAL(0,dp)
1473
      END DO
1474
      !$OMP END DO
1475
    END IF
1476
#endif
1477

1478
    !$OMP SINGLE
1479
    A % Rows(1) = 1
1480
    DO i=2,N
1481
       j = InvPerm((i-2)/Ndeg+1)
1482
       A % Rows(i) = A % Rows(i-1) + Ndeg*RowNonzeros(j)
1483
    END DO
1484
    j = InvPerm((n-1)/ndeg+1)
1485
    A % Rows(n+1) = A % Rows(N)  +  Ndeg*RowNonzeros(j)
1486
    !$OMP END SINGLE
1487
    
1488
#ifdef _OPENMP
1489
    ! First touch matrix values with similar access pattern as in sparse dgemv
1490
    !$OMP DO
1491
    DO i=1,A % NumberOfRows
1492
      DO j=A % Rows(i), A % Rows(i+1)-1
1493
        A % Cols(j) = REAL(0,dp)
1494
      END DO
1495
    END DO
1496
    !$OMP END DO NOWAIT
1497
    !$OMP DO
1498
    DO i=1,n
1499
      A % Diag(i) = REAL(0,dp)
1500
    END DO
1501
    !$OMP END DO
1502
#else
1503
    A % Cols = 0
1504
    A % Diag = 0
1505
#endif
1506
    !$OMP END PARALLEL
1507
    
1508

1509
    CALL Info('CRS_CreateMatrix','Creating CRS Matrix finished',Level=14)
1510

1511
  END FUNCTION CRS_CreateMatrix
1512
!------------------------------------------------------------------------------
1513

1514

1515
!------------------------------------------------------------------------------
1516
!>    Matrix vector product (v = Au) for a matrix given in CRS format.
1517
!------------------------------------------------------------------------------
1518
  SUBROUTINE CRS_MatrixVectorMultiply( A,u,v )
1519
!------------------------------------------------------------------------------
1520
    REAL(KIND=dp), DIMENSION(*), INTENT(IN) :: u   !< Vector to be multiplied
1521
    REAL(KIND=dp), DIMENSION(*), INTENT(OUT) :: v  !< Result vector
1522
    TYPE(Matrix_t), INTENT(IN) :: A                !< Structure holding matrix
1523
!------------------------------------------------------------------------------
1524
     INTEGER, POINTER  CONTIG :: Cols(:),Rows(:)
1525
     REAL(KIND=dp), POINTER  CONTIG :: Values(:)
1526

1527
     INTEGER :: i,j,n,k,l,m
1528
     REAL(KIND=dp) :: r1,r2,r3,r4,r5
1529
#ifdef HAVE_MKL
1530
     INTERFACE
1531
       SUBROUTINE mkl_dcsrgemv(transa, m, a, ia, ja, x, y)
1532
         USE Types
1533
         CHARACTER :: transa
1534
         INTEGER :: m
1535
         REAL(KIND=dp) :: a(*)
1536
         INTEGER :: ia(*), ja(*)
1537
         REAL(KIND=dp) :: x(*), y(*)
1538
       END SUBROUTINE mkl_dcsrgemv
1539
     END INTERFACE
1540
#endif
1541

1542
!------------------------------------------------------------------------------
1543

1544
    n = A % NumberOfRows
1545
    Rows   => A % Rows
1546
    Cols   => A % Cols
1547
    Values => A % Values
1548
    
1549
    IF  ( A % MatvecSubr /= 0 ) THEN
1550
      CALL MatVecSubrExt(A % MatVecSubr,A % SpMV, n,Rows,Cols,Values,u,v,0)
1551
      RETURN
1552
    END IF
1553

1554
    ! Use MKL to perform mvp if it is available
1555
#ifdef HAVE_MKL
1556
    CALL mkl_dcsrgemv('N', n, Values, Rows, Cols, u, v)
1557
#else
1558

1559
    ! There may be a small structured block in the CRS matrix that is due to the problem
1560
    ! being initially vector valued. For example, in 3D elasticity we usually have dofs related
1561
    ! to (x,y,z) displacements following each other. Using this small dense block we may reduce
1562
    ! indirect memory addressing a little.
1563
    !-------------------------------------------------------------------------------------------
1564
    SELECT CASE( A % ndeg )
1565
      
1566
    CASE( 5, 10 )
1567
      !$omp parallel do private(j,l,r1,r2,r3,r4,r5)
1568
      DO i=1,n
1569
        r1 = 0.0_dp; r2 = 0.0_dp; r3 = 0.0_dp; r4 = 0.0_dp; r5 = 0.0_dp
1570
        !DIR$ IVDEP
1571
        DO j=Rows(i),Rows(i+1)-1,5
1572
          l = Cols(j)
1573
          r1 = r1 + u(l) * Values(j)
1574
          r2 = r2 + u(l+1) * Values(j+1)
1575
          r3 = r3 + u(l+2) * Values(j+2)
1576
          r4 = r4 + u(l+3) * Values(j+3)
1577
          r5 = r5 + u(l+4) * Values(j+4)
1578
        END DO
1579
        v(i) = r1 + r2 + r3 + r4 + r5
1580
      END DO
1581
      !$omp end parallel do
1582

1583
    CASE( 4, 8 )
1584
      !$omp parallel do private(j,l,r1,r2,r3,r4)
1585
      DO i=1,n
1586
        r1 = 0.0_dp; r2 = 0.0_dp; r3 = 0.0_dp; r4 = 0.0_dp
1587
        !DIR$ IVDEP
1588
        DO j=Rows(i),Rows(i+1)-1,4
1589
          l = Cols(j)
1590
          r1 = r1 + u(l) * Values(j)
1591
          r2 = r2 + u(l+1) * Values(j+1)
1592
          r3 = r3 + u(l+2) * Values(j+2)
1593
          r4 = r4 + u(l+3) * Values(j+3)
1594
        END DO
1595
        v(i) = r1 + r2 + r3 + r4
1596
      END DO
1597
      !$omp end parallel do
1598

1599
    CASE( 3, 6 )
1600
      !$omp parallel do private(j,l,r1,r2,r3)
1601
      DO i=1,n
1602
        r1 = 0.0_dp; r2 = 0.0_dp; r3 = 0.0_dp
1603
        !DIR$ IVDEP
1604
        DO j=Rows(i),Rows(i+1)-1,3
1605
          l = Cols(j)
1606
          r1 = r1 + u(l) * Values(j)
1607
          r2 = r2 + u(l+1) * Values(j+1)
1608
          r3 = r3 + u(l+2) * Values(j+2)
1609
        END DO
1610
        v(i) = r1 + r2 + r3
1611
      END DO
1612
      !$omp end parallel do
1613
      
1614
    CASE( 2 )
1615
      !$omp parallel do private(j,l,r1,r2)
1616
      DO i=1,n
1617
        r1 = 0.0_dp; r2 = 0.0_dp
1618
        !DIR$ IVDEP
1619
        DO j=Rows(i),Rows(i+1)-1,2
1620
          l = Cols(j)
1621
          r1 = r1 + u(l) * Values(j)
1622
          r2 = r2 + u(l+1) * Values(j+1)
1623
        END DO
1624
        v(i) = r1 + r2
1625
      END DO
1626
      !$omp end parallel do
1627
      
1628
    CASE DEFAULT      
1629
      !$omp parallel do private(j,r1)
1630
      DO i=1,n
1631
        r1 = 0.0_dp
1632
        !DIR$ IVDEP
1633
        DO j=Rows(i),Rows(i+1)-1
1634
          r1 = r1 + u(Cols(j)) * Values(j)
1635
        END DO
1636
        v(i) = r1 
1637
      END DO
1638
      !$omp end parallel do      
1639
    END SELECT
1640
#endif
1641
!------------------------------------------------------------------------------
1642
  END SUBROUTINE CRS_MatrixVectorMultiply
1643
!------------------------------------------------------------------------------
1644

1645

1646
!------------------------------------------------------------------------------
1647
!> Matrix-vector multiply without initializing v to zero.
1648
!------------------------------------------------------------------------------
1649
  SUBROUTINE CRS_AdditiveMatrixVectorMultiply( A,u,v,c )
1650
!------------------------------------------------------------------------------
1651
    REAL(KIND=dp), DIMENSION(*), INTENT(IN) :: u   !< Vector to be multiplied
1652
    REAL(KIND=dp), DIMENSION(*), INTENT(OUT) :: v  !< Result vector
1653
    TYPE(Matrix_t), INTENT(IN) :: A                !< Structure holding matrix
1654
    REAL(KIND=dp), OPTIONAL :: c                   !< multiplier    
1655
    !------------------------------------------------------------------------------
1656
     INTEGER, POINTER  CONTIG :: Cols(:),Rows(:)
1657
     REAL(KIND=dp), POINTER  CONTIG :: Values(:)
1658
     INTEGER :: i,j,n
1659
     REAL(KIND=dp) :: rsum
1660

1661
!------------------------------------------------------------------------------
1662

1663
     n = A % NumberOfRows
1664
     Rows   => A % Rows
1665
     Cols   => A % Cols
1666
     Values => A % Values
1667

1668
     !$omp parallel do private(j,rsum)
1669
     DO i=1,n
1670
       rsum = 0.0d0
1671
       !DIR$ IVDEP
1672
       DO j=Rows(i),Rows(i+1)-1
1673
         rsum = rsum + u(Cols(j)) * Values(j)
1674
       END DO
1675
       
1676
       IF( PRESENT(c) ) THEN
1677
         v(i) = v(i) + c * rsum
1678
       ELSE
1679
         v(i) = v(i) + rsum
1680
       END IF
1681
     END DO
1682
     !$omp end parallel do
1683
!------------------------------------------------------------------------------
1684
   END SUBROUTINE CRS_AdditiveMatrixVectorMultiply
1685
!------------------------------------------------------------------------------
1686

1687
   
1688
!------------------------------------------------------------------------------
1689
!> Matrix vector product (v = Au) for a matrix given in CRS format
1690
!> This one only applies to the active elements of u. The idea is that
1691
!> we may look at the partial matrix norm, for example. 
1692
!------------------------------------------------------------------------------
1693
  SUBROUTINE CRS_MaskedMatrixVectorMultiply( A,u,v,ActiveRow, ActiveCol )
1694
!------------------------------------------------------------------------------
1695
    REAL(KIND=dp), DIMENSION(*), INTENT(IN) :: u   !< Vector to be multiplied
1696
    REAL(KIND=dp), DIMENSION(*), INTENT(OUT) :: v  !< Result vector
1697
    TYPE(Matrix_t), INTENT(IN) :: A                !< Structure holding matrix
1698
    LOGICAL, DIMENSION(*), INTENT(IN) :: ActiveRow(:) !< Vector giving the active rows
1699
    LOGICAL, DIMENSION(*), INTENT(IN) :: ActiveCol(:) !< Vector giving the active columns
1700
    !------------------------------------------------------------------------------
1701
    INTEGER, POINTER  CONTIG :: Cols(:),Rows(:)
1702
    REAL(KIND=dp), POINTER  CONTIG :: Values(:)
1703
    INTEGER :: i,j,k,n
1704
    REAL(KIND=dp) :: rsum
1705

1706
    !------------------------------------------------------------------------------
1707

1708
    n = A % NumberOfRows
1709
    Rows   => A % Rows
1710
    Cols   => A % Cols
1711
    Values => A % Values
1712

1713
    DO i=1,n
1714
      IF( ActiveRow(i) ) THEN
1715
        rsum = 0.0d0
1716
        DO j=Rows(i),Rows(i+1)-1
1717
          k = Cols(j)
1718
          IF( ActiveCol(k) ) THEN
1719
            rsum = rsum + u(k) * Values(j)
1720
          END IF
1721
        END DO
1722
        v(i) = rsum
1723
      ELSE
1724
        v(i) = 0.0_dp
1725
      END IF
1726
    END DO
1727
!------------------------------------------------------------------------------
1728
  END SUBROUTINE CRS_MaskedMatrixVectorMultiply
1729
!------------------------------------------------------------------------------
1730

1731

1732
!------------------------------------------------------------------------------
1733
!> Matrix vector product (v = Au) for a matrix given in CRS format
1734
!> This one only applies to the active elements of u. The idea is that
1735
!> we may look at the partial matrix norm, for example. 
1736
!------------------------------------------------------------------------------
1737
  FUNCTION CRS_MatrixRowVectorMultiply( A,u,i) RESULT ( rsum )
1738
!------------------------------------------------------------------------------
1739
    TYPE(Matrix_t), INTENT(IN) :: A                !< Structure holding matrix
1740
    REAL(KIND=dp), DIMENSION(*), INTENT(IN) :: u   !< Vector to be multiplied
1741
    INTEGER :: i                                   !< Row of matrix
1742
    REAL(KIND=dp) :: rsum                          !< Matrix row x vector sum. 
1743
    !------------------------------------------------------------------------------
1744
    INTEGER, POINTER  CONTIG :: Cols(:),Rows(:)
1745
    REAL(KIND=dp), POINTER  CONTIG :: Values(:)
1746
    INTEGER :: j,n
1747
    !------------------------------------------------------------------------------
1748

1749
    IF(i<1 .OR. i>A % NumberOfRows) THEN
1750
      CALL Fatal('CRS_MatrixRowVectorMultiply','Invalid row number '//I2S(i))
1751
    END IF
1752
    Rows   => A % Rows
1753
    Cols   => A % Cols
1754
    Values => A % Values
1755

1756
    rsum = 0.0_dp
1757
    DO j=Rows(i),Rows(i+1)-1
1758
      rsum = rsum + Values(j) * u(Cols(j))
1759
    END DO 
1760
!------------------------------------------------------------------------------
1761
  END FUNCTION CRS_MatrixRowVectorMultiply
1762
!------------------------------------------------------------------------------
1763

1764
  
1765
  
1766
!------------------------------------------------------------------------------
1767
!>  Matrix-vector product v = |A|u with A a matrix in the CRS format and
1768
!>  |.| the matrix function giving the absolute values of the argument 
1769
!>  components. This special mv subroutine may be needed in connection with
1770
!>  certain stopping criteria for iterative linear solvers. 
1771
!------------------------------------------------------------------------------
1772
  SUBROUTINE CRS_ABSMatrixVectorMultiply( A,u,v )
1773
!------------------------------------------------------------------------------
1774
    REAL(KIND=dp), DIMENSION(*), INTENT(IN) :: u   !< The vector u
1775
    REAL(KIND=dp), DIMENSION(*), INTENT(OUT) :: v  !< The result vector v
1776
    TYPE(Matrix_t), INTENT(IN) :: A                !< The structure holding the matrix A
1777
!------------------------------------------------------------------------------
1778
    INTEGER, POINTER  CONTIG :: Cols(:),Rows(:)
1779
    REAL(KIND=dp), POINTER  CONTIG :: Values(:), Abs_Values(:)
1780

1781
    
1782
    INTEGER :: i,j,n
1783
    REAL(KIND=dp) :: rsum
1784
!------------------------------------------------------------------------------
1785

1786
    n = A % NumberOfRows
1787
    Rows   => A % Rows
1788
    Cols   => A % Cols
1789
    Values => A % Values
1790

1791
    IF  ( A % MatvecSubr /= 0 ) THEN
1792
      ALLOCATE(Abs_Values(SIZE(A % Values)))
1793
      Abs_Values = ABS(Values)
1794
      CALL MatVecSubrExt(A % MatVecSubr,A % SpMV, n,Rows,Cols,Abs_Values,u,v,0) ! TODO: (bug) must be ABS(Values)
1795
      DEALLOCATE(Abs_Values)
1796
      RETURN
1797
    END IF
1798

1799
!$omp parallel do private(j,rsum)
1800
    DO i=1,n
1801
      rsum = 0.0d0
1802
!DIR$ IVDEP
1803
      DO j=Rows(i),Rows(i+1)-1
1804
        rsum = rsum + u(Cols(j)) * ABS(Values(j))
1805
      END DO
1806
      v(i) = rsum
1807
    END DO
1808
!$omp end parallel do
1809
!------------------------------------------------------------------------------
1810
  END SUBROUTINE CRS_ABSMatrixVectorMultiply
1811
!------------------------------------------------------------------------------
1812

1813

1814
!------------------------------------------------------------------------------
1815
!>  Calculate transpose of A in CRS format: B = A^T
1816
!------------------------------------------------------------------------------
1817
     FUNCTION CRS_Transpose( A ) RESULT(B)
1818
!------------------------------------------------------------------------------
1819
       IMPLICIT NONE
1820
       
1821
       TYPE(Matrix_t), POINTER :: A, B
1822
       
1823
       INTEGER, ALLOCATABLE :: Row(:)
1824
       INTEGER :: NVals
1825
       INTEGER :: i,j,k,istat,nb, na
1826

1827
       CALL Info('CRS_Transpose','Creating a transpose of matrix',Level=20)
1828
       
1829
       B => AllocateMatrix()
1830

1831
       IF(.NOT. ASSOCIATED(A) ) THEN
1832
         CALL Fatal('CRS_Transpose','Matrix not associated!')
1833
       END IF
1834
         
1835
       na = A % NumberOfRows
1836
       IF( na == 0 ) THEN
1837
         B % NumberOfRows = 0
1838
         RETURN
1839
       END IF
1840
       
1841
       NVals = SIZE( A % Values )
1842
       nb = MAXVAL( A % Cols )
1843
       B % NumberOfRows = nb
1844
       
1845
       ALLOCATE( B % Rows( nb +1 ), B % Cols( NVals ), &
1846
           B % Values( Nvals ), Row( nb ), STAT=istat )
1847
       IF ( istat /= 0 )  CALL Fatal( 'CRS_Transpose','Memory allocation error.' )
1848

1849
       IF( ASSOCIATED( A % Diag ) ) THEN
1850
         ALLOCATE( B % Diag(nb) )       
1851
         B % Diag = 0
1852
       END IF
1853

1854
       ! Count how many hits there are in A for each column
1855
       Row = 0       
1856
       DO i = 1, NVals
1857
         Row( A % Cols(i) ) = Row( A % Cols(i) ) + 1
1858
       END DO
1859

1860
       ! For transpose the column hits are row hits
1861
       B % Rows = 0
1862
       B % Rows(1) = 1
1863
       DO i = 1, nB
1864
         B % Rows(i+1) = B % Rows(i) + Row(i)
1865
       END DO
1866
       B % Cols = 0
1867
       
1868
       ! Location of 1st entry in each row
1869
       Row(1:nB) = B % Rows(1:nB)
1870
       
1871
       DO i = 1, nA
1872

1873
         DO j = A % Rows(i), A % Rows(i+1) - 1
1874
           k = A % Cols(j)
1875

1876
           IF ( Row(k) < B % Rows(k+1) ) THEN 
1877
             B % Cols( Row(k) ) = i
1878
             B % Values( Row(k) ) = A % Values(j)
1879
             Row(k) = Row(k) + 1
1880
           ELSE
1881
             WRITE( Message, * ) 'Trying to access column beyond allocation: ', i,k,j
1882
             CALL Error( 'CRS_Transpose', Message )
1883
             RETURN
1884
           END IF
1885
         END DO
1886
       END DO
1887
              
1888
       DEALLOCATE( Row )
1889
       
1890
!------------------------------------------------------------------------------
1891
     END FUNCTION CRS_Transpose
1892
!------------------------------------------------------------------------------
1893

1894
!------------------------------------------------------------------------------
1895
!>    Matrix vector product (v = A^T u) for a matrix given in CRS format.
1896
!------------------------------------------------------------------------------
1897
  SUBROUTINE CRS_TransposeMatrixVectorMultiply( A,u,v )
1898
!------------------------------------------------------------------------------
1899
    REAL(KIND=dp), DIMENSION(*), INTENT(IN) :: u   !< Vector to be multiplied
1900
    REAL(KIND=dp), DIMENSION(*), INTENT(OUT) :: v  !< Result vector
1901
    TYPE(Matrix_t), INTENT(IN) :: A                !< Structure holding matrix
1902
!------------------------------------------------------------------------------
1903
     INTEGER, POINTER  CONTIG :: Cols(:),Rows(:)
1904
     REAL(KIND=dp), POINTER  CONTIG :: Values(:)
1905

1906
     INTEGER :: i,j,k,n,m
1907
     REAL(KIND=dp) :: rsum
1908
#ifdef HAVE_MKL
1909
	INTERFACE
1910
		SUBROUTINE mkl_dcsrgemv(transa, m, a, ia, ja, x, y)
1911
	 		USE Types
1912
	 		CHARACTER :: transa
1913
	 		INTEGER :: m
1914
	 		REAL(KIND=dp) :: a(*)
1915
	 		INTEGER :: ia(*), ja(*)
1916
	 		REAL(KIND=dp) :: x(*), y(*)
1917
	 	END SUBROUTINE mkl_dcsrgemv
1918
	END INTERFACE
1919
#endif
1920
!------------------------------------------------------------------------------
1921

1922
     n = A % NumberOfRows
1923
     Rows   => A % Rows
1924
     Cols   => A % Cols
1925
     Values => A % Values     
1926
          
1927
     ! Use MKL to perform mvp if it is available
1928
#ifdef HAVE_MKL
1929
     CALL mkl_dcsrgemv('T', n, Values, Rows, Cols, u, v)
1930
#else
1931
     m = MAXVAL(Cols)
1932
     v(1:m) = 0.0_dp
1933
     DO i=1,n
1934
!DIR$ IVDEP
1935
       DO j=Rows(i),Rows(i+1)-1
1936
         k = Cols(j)
1937
         v(k) = v(k) + u(i) * Values(j)
1938
       END DO
1939
     END DO
1940
#endif
1941
!------------------------------------------------------------------------------
1942
   END SUBROUTINE CRS_TransposeMatrixVectorMultiply
1943
!------------------------------------------------------------------------------
1944

1945

1946
!------------------------------------------------------------------------------
1947
!> Add another matrix B to matrix A, or created a combined matrix C. 
1948
!------------------------------------------------------------------------------
1949
  SUBROUTINE CRS_MergeMatrix( A,B,C,PermA,PermB,PermC)
1950
!------------------------------------------------------------------------------
1951
    TYPE(Matrix_t), POINTER :: A             !< Structure holding the master matrix
1952
    TYPE(Matrix_t), POINTER :: B             !< Structure holding the slave matrix
1953
    TYPE(Matrix_t), POINTER, OPTIONAL :: C   !< Structure holding the sum matrix
1954
    INTEGER, POINTER, OPTIONAL :: PermA(:)   !< Permutation of the master dofs
1955
    INTEGER, POINTER, OPTIONAL :: PermB(:)   !< Permutation of the slave dofs
1956
    INTEGER, POINTER, OPTIONAL :: PermC(:)   !< Permutation of the combined dofs
1957
!------------------------------------------------------------------------------
1958
    INTEGER, POINTER  CONTIG :: ColsA(:),RowsA(:),ColsB(:),RowsB(:),&
1959
        Rows(:),Cols(:),invPermA(:),invPermB(:),invPermC(:),Perm(:)
1960
    REAL(KIND=dp), POINTER  CONTIG :: ValuesA(:),ValuesB(:),Values(:),rhs(:)
1961
    INTEGER :: i,j,k,n,m,nA,nB,nC,kb,kb0,iA,iB,iC,colj
1962
    LOGICAL :: Set,UsePerm
1963
    INTEGER, ALLOCATABLE :: ColUsed(:)
1964
!------------------------------------------------------------------------------
1965

1966
     CALL Info('CRS_MergeMatrix','Merging two matrices',Level=9)
1967

1968
     IF(.NOT. ASSOCIATED(A)) THEN
1969
       CALL Fatal('CRS_MergeMatrix','A not associated')
1970
     ELSE IF(.NOT. ASSOCIATED(B)) THEN
1971
       CALL Fatal('CRS_MergeMatrix','B not associated')
1972
     END IF
1973

1974
     UsePerm = PRESENT( PermA ) 
1975
     
1976
     IF( UsePerm ) THEN
1977
       IF(.NOT. PRESENT( PermB ) ) THEN
1978
         CALL Fatal('CRS_MergeMatrix','Either both PermA and PermB or neither')
1979
       END IF        
1980
       n = SIZE(PermA)
1981
       IF( SIZE(PermB) /= n ) THEN
1982
         CALL Fatal('CRS_MergeMatrix','Mismatch in perm size')
1983
       END IF
1984
     ELSE
1985
       n = MAX( A % NumberOfRows, B % NumberOfRows ) 
1986
     END IF
1987

1988
     RowsA   => A % Rows
1989
     ColsA   => A % Cols
1990
     ValuesA => A % Values
1991
            
1992
     RowsB   => B % Rows
1993
     ColsB   => B % Cols
1994
     ValuesB => B % Values
1995
     Set = .FALSE.
1996
     
1997
     IF( UsePerm ) THEN
1998
       nA = MAXVAL( permA )
1999
       nB = MAXVAL( permB )
2000

2001
       ALLOCATE(InvPermA(na),InvPermB(nB))
2002
       invPermA = 0
2003
       invPermB = 0
2004
       DO i=1,SIZE(permA)
2005
         IF( permA(i) > 0) invPermA(permA(i)) = i
2006
         IF( permB(i) > 0) invPermB(permB(i)) = i
2007
       END DO
2008

2009
       NULLIFY(Perm)
2010
       ALLOCATE(Perm(SIZE(permA)))
2011
       IF(PRESENT(PermC)) PermC => Perm
2012
       Perm = PermA       
2013
       j = MAXVAL(PermA)
2014

2015
       DO i=1,SIZE(PermB)
2016
         IF(PermA(i) == 0 .AND. PermB(i) > 0 ) THEN
2017
           j = j+1
2018
           Perm(i) = j
2019
         END IF
2020
       END DO
2021

2022
       ! Fast way to check whether the entry has been created.
2023
       ALLOCATE(ColUsed(j)) 
2024
       ColUsed = 0
2025
     END IF
2026
       
2027
100  kb = 0
2028
     iC = 0
2029
     IF( UsePerm ) THEN         
2030
       DO iC=1,SIZE(InvPermA)
2031
         i = InvPermA(iC)
2032

2033
         iA = PermA(i)
2034
         IF(iA /= iC) CALL Fatal('CRS_MergeMatrix','This Should be True by consruction!')                 
2035
         iB = PermB(i)
2036
         
2037
         nA = 0
2038
         IF( iA > 0 ) nA = RowsA(iA+1)-RowsA(iA) 
2039
         nB = 0
2040
         IF( iB > 0 ) nB = RowsB(iB+1)-RowsB(iB)         
2041

2042
         IF(nA == 0) CALL Fatal('CRS_MergeMatrix','This should not happen for nA!')
2043

2044
         ! Do the case with A active (with or without B)
2045
         IF( nB > 0 ) THEN
2046
           DO j=RowsA(iA),RowsA(iA+1)-1
2047
             kb = kb + 1
2048
             colj = Perm(invPermA(ColsA(j)))
2049
             ColUsed(colj) = kb
2050

2051
             IF( Set ) THEN
2052
               Cols(kb) = colj
2053
               Values(kb) = ValuesA(j)
2054
             END IF
2055
           END DO           
2056
           DO j=RowsB(iB),RowsB(iB+1)-1
2057
             colj = Perm(invPermB(ColsB(j)))
2058
             kb0 = ColUsed(colj)
2059

2060
             IF(kb0 > 0 ) THEN
2061
               IF( Set ) THEN
2062
                 Values(kb0) = Values(kb0) + ValuesB(j)
2063
               END IF
2064
             ELSE
2065
               kb = kb + 1
2066
               IF( Set ) THEN
2067
                 Cols(kb) = colj
2068
                 Values(kb) = ValuesB(j)
2069
               END IF
2070
             END IF
2071
           END DO
2072
           DO j=RowsA(iA),RowsA(iA+1)-1
2073
             colj = Perm(invPermA(ColsA(j)))
2074
             ColUsed(colj) = 0
2075
           END DO
2076
         ELSE
2077
           DO j=RowsA(iA),RowsA(iA+1)-1
2078
             kb = kb + 1
2079
             IF( Set ) THEN
2080
               colj = Perm(invPermA(ColsA(j)))
2081
               Cols(kb) = colj
2082
               Values(kb) = ValuesA(j)
2083
             END IF
2084
           END DO
2085
         END IF
2086
         IF( Set ) THEN
2087
           Rows(iC+1) = kb+1
2088
         END IF
2089
       END DO
2090
       
2091
       ! Do the nodes with only B active
2092
       iC = SIZE(InvPermA)
2093
       DO i=1,n
2094
         iA = PermA(i)
2095
         iB = PermB(i)
2096

2097
         IF(iA > 0 .OR. iB == 0) CYCLE
2098

2099
         nB = RowsB(iB+1)-RowsB(iB)         
2100
         iC = Perm(i)
2101

2102
         DO j=RowsB(iB),RowsB(iB+1)-1
2103
           kb = kb + 1
2104
           IF( Set ) THEN
2105
             colj = Perm(invPermB(ColsB(j)))
2106
             Cols(kb) = colj
2107
             Values(kb) = ValuesB(j)
2108
           END IF
2109
         END DO
2110
         IF( Set ) THEN
2111
           Rows(iC+1) = kb+1
2112
         END IF
2113
       END DO
2114
     ELSE
2115
       DO i=1,n
2116
         nA = RowsA(i+1)-RowsA(i) 
2117
         nB = RowsB(i+1)-RowsB(i)
2118
         IF( nA > 0 .AND. nB > 0 ) THEN
2119
           WRITE (Message,'(A,I0,I0)') 'Code the possibility to merge rows: ',iA,iB
2120
           CALL Fatal('CRS_MergeMatrix',Message)
2121
         ELSE IF( nA > 0 ) THEN
2122
           DO j=RowsA(i),RowsA(i+1)-1
2123
             kb = kb + 1
2124
             IF( Set ) THEN
2125
               Cols(kb) = ColsA(j)
2126
               Values(kb) = ValuesA(j)
2127
             END IF
2128
           END DO
2129
         ELSE IF( nB > 0 ) THEN
2130
           DO j=RowsB(i),RowsB(i+1)-1
2131
             kb = kb + 1
2132
             IF( Set ) THEN
2133
               Cols(kb) = ColsB(j)
2134
               Values(kb) = ValuesB(j)
2135
             END IF
2136
           END DO
2137
         END IF
2138
         IF( Set ) THEN
2139
           Rows(i+1) = kb+1
2140
         END IF
2141
       END DO
2142
     END IF
2143

2144
     IF( kb == 0 ) THEN
2145
       CALL Fatal('CRS_MergeMatrix','Union size is zero?')
2146
     END IF
2147
     
2148
     IF(.NOT. Set) THEN
2149
       IF( UsePerm ) THEN
2150
         nC = iC
2151
       ELSE
2152
         nC = n
2153
       END IF       
2154
       ALLOCATE( Rows(nC+1), Cols(kb), Values(kb) )
2155
       Rows = 0
2156
       Cols = 0
2157
       Values = 0.0_dp
2158
       Rows(1) = 1
2159

2160
       CALL Info('CRS_MergeMatrix','Combined matrix has '//I2S(nC)//' rows',Level=10)
2161
       CALL Info('CRS_MergeMatrix','Combined matrix has '//I2S(kb)//' nonzeros',Level=10)
2162
       
2163
       Set = .TRUE.
2164
       CALL Info('CRS_MergeMatrix','Done Allocating and going now really',Level=9)
2165
       GOTO 100
2166
     END IF
2167
     
2168
     IF( PRESENT(C) ) THEN
2169
       C % Rows => Rows
2170
       C % Cols => Cols
2171
       C % Values => Values
2172
       C % NumberOfRows = nC
2173
     ELSE              
2174
       DEALLOCATE( RowsA, RowsB, ColsA, ColsB, ValuesA, ValuesB )     
2175
       B % NumberOfRows = 0       
2176
       A % Rows => Rows
2177
       A % Cols => Cols
2178
       A % Values => Values
2179
       A % NumberOfRows = nC
2180
     END IF
2181

2182
     IF(UsePerm) THEN
2183
       DEALLOCATE(invPermA, invPermB, ColUsed)
2184
     END IF
2185
    
2186
     CALL Info('CRS_MergeMatrix','Merging of matrices finished',Level=9)
2187

2188
!------------------------------------------------------------------------------
2189
   END SUBROUTINE CRS_MergeMatrix
2190
!------------------------------------------------------------------------------
2191

2192

2193
!------------------------------------------------------------------------------
2194
!>    Matrix vector product (v = Au) for a matrix given in CRS format
2195
!>    assuming complex valued matrix equation.
2196
!------------------------------------------------------------------------------
2197
  SUBROUTINE CRS_ComplexMatrixVectorMultiply( A,u,v )
2198
!------------------------------------------------------------------------------
2199
    COMPLEX(KIND=dp), DIMENSION(*), INTENT(IN) :: u   !< Vector to be multiplied
2200
    COMPLEX(KIND=dp), DIMENSION(*), INTENT(OUT) :: v  !< Result vector
2201
    TYPE(Matrix_t), INTENT(IN) :: A                !< Structure holding matrix
2202
!------------------------------------------------------------------------------
2203
    INTEGER, POINTER :: Cols(:),Rows(:)
2204
    REAL(KIND=dp), POINTER :: Values(:)
2205
    INTEGER :: i,j,n
2206
    COMPLEX(KIND=dp) :: s,rsum
2207
!------------------------------------------------------------------------------
2208
    n = A % NumberOfRows / 2
2209
    Rows   => A % Rows
2210
    Cols   => A % Cols
2211
    Values => A % Values
2212

2213
!$omp parallel do private(rsum,j,s)
2214
    DO i=1,n
2215
       rsum = CMPLX( 0.0d0, 0.0d0,KIND=dp )
2216
!DIR$ IVDEP
2217
       DO j=Rows(2*i-1),Rows(2*i)-1,2
2218
          s = CMPLX( Values(j), -Values(j+1), KIND=dp )
2219
          rsum = rsum + s * u((Cols(j)+1)/2)
2220
       END DO
2221
       v(i) = rsum
2222
    END DO
2223
!$omp end parallel do
2224
!------------------------------------------------------------------------------
2225
  END SUBROUTINE CRS_ComplexMatrixVectorMultiply
2226
!------------------------------------------------------------------------------
2227

2228

2229
!-------------------------------------------------------------------------------
2230
  SUBROUTINE CRS_ApplyProjector( PMatrix, u, uperm, v, vperm, Trans )
2231
!-------------------------------------------------------------------------------
2232
    TYPE(Matrix_t) :: PMatrix
2233
    REAL(KIND=dp) :: u(:),v(:)
2234
    LOGICAL, OPTIONAL :: Trans
2235
    INTEGER, POINTER :: uperm(:), vperm(:)
2236
!-------------------------------------------------------------------------------
2237
    INTEGER :: i,j,k,l,n
2238
    REAL(KIND=dp), POINTER :: Values(:)
2239
    LOGICAL :: LTrans
2240
    INTEGER, POINTER :: Rows(:), Cols(:)
2241
!-------------------------------------------------------------------------------
2242
    LTrans = .FALSE.
2243
    IF ( PRESENT( Trans ) ) LTrans = Trans
2244

2245
    n = PMatrix % NumberOfRows
2246
    Rows   => PMatrix % Rows
2247
    Cols   => PMatrix % Cols
2248
    Values => PMatrix % Values
2249

2250
    IF ( ASSOCIATED( uperm ) .AND. ASSOCIATED( vperm ) ) THEN
2251
       IF ( LTrans ) THEN
2252
          DO i=1,n
2253
             k = uperm(i)
2254
             IF ( k > 0 ) THEN
2255
                DO j=Rows(i),Rows(i+1)-1
2256
                   l = vperm(Cols(j))
2257
                   IF ( l > 0 ) v(l) = v(l) + u(k) * Values(j)
2258
                END DO
2259
             END IF
2260
          END DO
2261
       ELSE
2262
          DO i=1,n
2263
             l = vperm(i)
2264
             IF ( l > 0 ) THEN
2265
               IF ( ANY(Values(Rows(i):Rows(i+1)-1)/=0) ) v(l)=0
2266
             END IF
2267
          END DO
2268

2269
          DO i=1,n
2270
             l = vperm(i)
2271
             IF ( l > 0 ) THEN
2272
                DO j = Rows(i), Rows(i+1)-1
2273
                   k = uperm(Cols(j))
2274
                   IF ( k>0 ) v(l) = v(l) + u(k) * Values(j)
2275
                END DO
2276
             END IF
2277
          END DO
2278
       END IF
2279
    ELSE
2280
       IF ( LTrans ) THEN
2281
          DO i=1,n
2282
             DO j=Rows(i),Rows(i+1)-1
2283
                v(Cols(j)) = v(Cols(j)) + u(i) * Values(j)
2284
             END DO
2285
          END DO
2286
       ELSE
2287
          DO i=1,n
2288
             DO j = Rows(i), Rows(i+1)-1
2289
                v(i) = v(i) + u(Cols(j)) * Values(j)
2290
             END DO
2291
          END DO
2292
       END IF
2293
    END IF
2294
!-------------------------------------------------------------------------------
2295
  END SUBROUTINE CRS_ApplyProjector
2296
!-------------------------------------------------------------------------------
2297

2298

2299

2300

2301
!------------------------------------------------------------------------------
2302
!>    Diagonal preconditioning of a CRS format matrix. Matrix is accessed
2303
!>    from a global variable GlobalMatrix. Note that if the matrix has been
2304
!> scaled so that the diagonal entries are already ones, this subroutine is obsolete.
2305
!------------------------------------------------------------------------------
2306
  SUBROUTINE CRS_DiagPrecondition( u,v,ipar )
2307
!------------------------------------------------------------------------------
2308
    REAL(KIND=dp), DIMENSION(*), INTENT(OUT) :: u  !< Resulting approximate solution after preconditioning
2309
    REAL(KIND=dp), DIMENSION(*), INTENT(IN) :: v  !< Given right-hand-side
2310
    INTEGER, DIMENSION(*) :: ipar   !< structure holding info from (HUTIter-iterative solver package)
2311
!------------------------------------------------------------------------------
2312
    INTEGER :: i,j,n
2313
    INTEGER, POINTER :: Cols(:),Rows(:),Diag(:)
2314
    REAL(KIND=dp), POINTER :: Values(:)
2315

2316
    Diag   => GlobalMatrix % Diag
2317
    Rows   => GlobalMatrix % Rows
2318
    Cols   => GlobalMatrix % Cols
2319
    Values => GlobalMatrix % Values
2320

2321
    n = GlobalMatrix % NumberOfRows
2322

2323
    IF ( .NOT. GlobalMatrix % Ordered ) THEN
2324
       !$OMP PARALLEL DO
2325
       DO i=1,N
2326
          CALL SortF( Rows(i+1)-Rows(i),Cols(Rows(i):Rows(i+1)-1), &
2327
                   Values(Rows(i):Rows(i+1)-1) )
2328
       END DO
2329
       !$OMP END PARALLEL DO
2330
       !$OMP PARALLEL DO
2331
       DO i=1,N
2332
          DO j=Rows(i),Rows(i+1)-1
2333
             IF ( Cols(j) == i ) THEN
2334
                Diag(i) = j
2335
                EXIT
2336
             END IF
2337
          END DO
2338
       END DO
2339
       !$OMP END PARALLEL DO
2340
       GlobalMatrix % Ordered = .TRUE.
2341
    END IF
2342

2343
    !$OMP PARALLEL DO
2344
    DO i=1,n
2345
       IF  ( ABS( Values(Diag(i))) > AEPS ) THEN
2346
           u(i) = v(i) / Values(Diag(i))
2347
       ELSE
2348
           u(i) = v(i)
2349
       END IF
2350
    END DO
2351
    !$OMP END PARALLEL DO
2352
!------------------------------------------------------------------------------
2353
  END SUBROUTINE CRS_DiagPrecondition
2354
!------------------------------------------------------------------------------
2355

2356

2357

2358
!------------------------------------------------------------------------------
2359
!>    Diagonal preconditioning of a CRS format matrix for complex valued matrix equations. 
2360
!>    Matrix is accessed from a global variable GlobalMatrix.
2361
!------------------------------------------------------------------------------
2362
  SUBROUTINE CRS_ComplexDiagPrecondition( u,v,ipar )
2363
!------------------------------------------------------------------------------
2364
    COMPLEX(KIND=dp), DIMENSION(*), INTENT(OUT) :: u  !< Resulting approximate solution after preconditioning
2365
    COMPLEX(KIND=dp), DIMENSION(*), INTENT(IN) :: v  !< Given right-hand-side
2366
    INTEGER, DIMENSION(*) :: ipar   !< structure holding info from (HUTIter-iterative solver package)
2367
!------------------------------------------------------------------------------
2368
    INTEGER :: i,j,n
2369
    INTEGER, POINTER :: Cols(:),Rows(:),Diag(:)
2370
    COMPLEX(KIND=dp) :: A
2371
    REAL(KIND=dp), POINTER :: Values(:)
2372

2373
    Diag   => GlobalMatrix % Diag
2374
    Rows   => GlobalMatrix % Rows
2375
    Cols   => GlobalMatrix % Cols
2376
    Values => GlobalMatrix % Values
2377

2378
    n = GlobalMatrix % NumberOfRows
2379

2380
    IF ( .NOT. GlobalMatrix % Ordered ) THEN
2381
       DO i=1,N
2382
          CALL SortF( Rows(i+1)-Rows(i),Cols(Rows(i):Rows(i+1)-1), &
2383
                   Values(Rows(i):Rows(i+1)-1) )
2384
       END DO
2385

2386
       DO i=1,N
2387
          DO j=Rows(i),Rows(i+1)-1
2388
             IF ( Cols(j) == i ) THEN
2389
                Diag(i) = j
2390
                EXIT
2391
             END IF
2392
          END DO
2393
       END DO
2394
       GlobalMatrix % Ordered = .TRUE.
2395
    END IF
2396

2397
    DO i=1,n/2
2398
       A = CMPLX( Values(Diag(2*i-1)), -Values(Diag(2*i-1)+1), KIND=dp )
2399
       u(i) = v(i) / A
2400
    END DO
2401
!------------------------------------------------------------------------------
2402
  END SUBROUTINE CRS_ComplexDiagPrecondition
2403
!------------------------------------------------------------------------------
2404

2405

2406
!------------------------------------------------------------------------------
2407
!> Picks the block diagonal entries from matrix A to build matrix B.
2408
!------------------------------------------------------------------------------
2409
  SUBROUTINE CRS_BlockDiagonal(A,B,Blocks) 
2410
!------------------------------------------------------------------------------
2411
    TYPE(Matrix_t), INTENT(IN) :: A      !< The initial matrix
2412
    TYPE(Matrix_t) :: B  !< The block diagonal matrix
2413
    INTEGER, INTENT(IN) :: Blocks        !< Number of blocks used in the decomposition
2414
!------------------------------------------------------------------------------
2415
    INTEGER :: i,j,k,l,kb,n
2416

2417
    IF(Blocks <= 1) RETURN
2418

2419
    N = A % NumberOfRows
2420
    B % NumberOfRows = N
2421
    
2422
    kb = 0
2423
    DO i=1,N
2424
      DO k= A % Rows(i), A % Rows(i+1)-1
2425
        l = A % Cols(k)
2426
        IF( MOD(i,Blocks) == MOD(l,Blocks)) kb = kb + 1
2427
      END DO
2428
    END DO
2429
    ALLOCATE(B % Rows(N+1),B % Cols(kb), B % Values(kb), B % Diag(n))
2430
      
2431
    kb = 1
2432
    DO i=1,N
2433
      B % Rows(i) = kb
2434
      DO k = A % Rows(i), A % Rows(i+1)-1
2435
        l = A % Cols(k)
2436
        IF( MOD(i,Blocks) == MOD(l,Blocks) ) THEN
2437
          B % Values(kb) = A % Values(k)
2438
          B % Cols(kb) = A % Cols(k)
2439
          IF( B % Cols(kb) == i) B % Diag(i) = kb
2440
          kb = kb + 1  
2441
        END IF
2442
      END DO
2443
    END DO
2444
    B % Rows(N+1) = kb
2445

2446
!------------------------------------------------------------------------------
2447
  END SUBROUTINE CRS_BlockDiagonal
2448
!------------------------------------------------------------------------------
2449

2450

2451
!------------------------------------------------------------------------------
2452
!> Removes zeros from the matrix structure.
2453
!> This might be done in order to save memory, or to speed up the matrix 
2454
!> operations. One must be careful since the fact the an entry is zero
2455
!> does not always imply that it would be zero throughout the simulation.
2456
!------------------------------------------------------------------------------
2457
  SUBROUTINE CRS_RemoveZeros( A, NoDiag, RemoveEps )
2458
!-------------------------------------------------------------------------------------------
2459
    TYPE(Matrix_t) :: A          !< The matrix which will be returned with the non-zeros removed
2460
    LOGICAL, OPTIONAL :: NoDiag  !< Can we also loose the diag if it happens to be zero? 
2461
    REAL(KIND=dp), OPTIONAL :: RemoveEps
2462
    !-------------------------------------------------------------------------------------------
2463
    INTEGER :: i,j,k,l,iml,kb,kb0,n,rowkb
2464
    INTEGER, POINTER CONTIG :: Cols(:), Rows(:), Diag(:)
2465
    REAL(KIND=DP) :: val, imval, reps
2466
    REAL(KIND=DP), POINTER CONTIG :: Values(:)
2467
    LOGICAL :: IsComplex, Hit, ImHit, CheckDiag
2468
    
2469
    N = A % NumberOfRows
2470

2471
    IsComplex = A % Complex 
2472

2473
    IF( PRESENT( NoDiag ) ) THEN
2474
      CheckDiag = .NOT. NoDiag
2475
    ELSE
2476
      CheckDiag = .TRUE.
2477
    END IF
2478

2479
    IF( PRESENT( RemoveEps ) ) THEN
2480
      reps = RemoveEps
2481
    ELSE
2482
      reps = ( EPSILON( reps ) ) **2 
2483
    END IF
2484
      
2485
    ! Count the number of nonzeros
2486
    ! The diagonal entry is assumed always to exist.
2487
    kb = 0
2488

2489
    IF( IsComplex ) THEN
2490
      DO i=1,N
2491
        DO k= A % Rows(i), A % Rows(i+1)-1, 2
2492
          l = A % Cols(k)
2493
          val = A % Values(k)
2494
          Hit = ( ( CheckDiag .AND. i == l ) .OR. ABS( val ) > reps )
2495

2496
          iml = A % Cols(k+1)
2497
          imval = A % Values(k+1)
2498
          ImHit = ( ( CheckDiag .AND. i == iml ) .OR. ABS( imval ) > reps )
2499
          
2500
          IF( Hit .OR. ImHit ) kb = kb + 2
2501
        END DO
2502
      END DO
2503
    ELSE
2504
      DO i=1,N
2505
        DO k= A % Rows(i), A % Rows(i+1)-1
2506
          l = A % Cols(k)
2507
          val = A % Values(k)
2508
          Hit = ( ( CheckDiag .AND. i == l ) .OR. ABS( val ) > reps )
2509
          IF( Hit ) kb = kb + 1
2510
        END DO
2511
      END DO
2512
    END IF
2513
      
2514
    kb0 = SIZE( A % Values )
2515

2516

2517
    IF( kb == kb0 ) THEN
2518
      CALL Info('CRS_RemoveZeros','There are no zeros to remove',Level=6)
2519
      RETURN
2520
    END IF
2521

2522
    WRITE( Message,'(A,F8.3,A)') 'Fraction of zeros to remove: ',100.0*(1.0-1.0*kb/kb0),' %'
2523
    CALL Info('CRS_RemoveZeros', Message)
2524

2525
    ! These are new
2526
    ALLOCATE(Cols(kb), Values(kb))
2527
   
2528
    ! These are overwritten
2529
    Diag => A % Diag
2530
    Rows => A % Rows
2531

2532

2533
    kb = 0
2534
    Rows(1) = 1
2535

2536
    IF( IsComplex ) THEN
2537
      DO i=1,N
2538
        kb0 = kb+1
2539
        DO k = A % Rows(i), A % Rows(i+1)-1,2
2540
          l = A % Cols(k)
2541
          val = A % Values(k) 
2542
          Hit = ( ( CheckDiag .AND. i == l ) .OR. ABS( val ) > reps )
2543

2544
          iml = A % Cols(k+1)
2545
          imval = A % Values(k+1)
2546
          ImHit = ( ( CheckDiag .AND. i == iml ) .OR. ABS( imval ) > reps ) 
2547

2548
          IF( Hit .OR. ImHit ) THEN
2549
            kb = kb + 1
2550
            IF( CheckDiag .AND. i == l ) Diag(i) = kb
2551
            Values(kb) = val
2552
            Cols(kb) = l
2553
            
2554
            kb = kb + 1
2555
            IF( CheckDiag .AND. i == iml ) Diag(i) = kb
2556
            Values(kb) = imval
2557
            Cols(kb) = iml
2558
          END IF
2559
        END DO
2560
        Rows(i) = kb0
2561
      END DO
2562
      
2563
    ELSE      
2564
      DO i=1,N
2565
        kb0 = kb+1
2566
        DO k = A % Rows(i), A % Rows(i+1)-1
2567
          l = A % Cols(k)
2568
          val = A % Values(k) 
2569

2570
          Hit = ( i == l .OR. ABS( val ) > reps )
2571
          
2572
          IF( Hit ) THEN
2573
            kb = kb + 1
2574
            IF( CheckDiag .AND. i == l ) Diag(i) = kb
2575
            
2576
            ! Set the new entry to the matrix
2577
            Values(kb) = val
2578
            Cols(kb) = l
2579
          END IF
2580
            
2581
        END DO
2582
        Rows(i) = kb0
2583
      END DO
2584
    END IF
2585
    Rows(N+1) = kb+1
2586
    
2587
    
2588
    DEALLOCATE( A % Values, A % Cols ) 
2589
    A % Values => Values
2590
    A % Cols => Cols
2591

2592
    ! This can no longer have structured blocks
2593
    A % Ndeg = -1
2594
    
2595
    IF(.NOT. CheckDiag ) THEN
2596
      IF( ASSOCIATED( A % Diag ) ) DEALLOCATE( A % Diag )
2597
    END IF
2598
    
2599
!------------------------------------------------------------------------------
2600
  END SUBROUTINE CRS_RemoveZeros
2601
!------------------------------------------------------------------------------
2602
 
2603
!------------------------------------------------------------------------------
2604
!> Makes a algebraic lower order scheme assuming steady state advection-diffusion equation.
2605
!> This can be applied together with flux corrected transport (FCT) scheme.
2606
!> Also creates a lumped mass to MassValuesLumped and saves the original stiffness
2607
!> matrix values to BulkValues.
2608
!
2609
!> For more information see, for example, 
2610
!> Dmitri Kuzmin (2008): "Explicit and implicit FEM-FCT algorithms with flux linearization"
2611
!------------------------------------------------------------------------------
2612
  SUBROUTINE CRS_FCTLowOrder( A )
2613
    USE SparIterGlobals
2614
#if defined(ELMER_HAVE_MPI_MODULE)
2615
    USE mpi
2616
#endif
2617
!------------------------------------------------------------------------------
2618
    TYPE(Matrix_t) :: A           !< Initial higher order matrix
2619
!------------------------------------------------------------------------------
2620
    INTEGER :: i,j,k,k2,l,kb,n, proc0, proc1,p,q, gi,gj
2621
    REAL(KIND=dp) :: Aij,Aji,Aii,Dij,dFij,msum
2622
    REAL(KIND=dp), POINTER :: ML(:)
2623
    INTEGER, POINTER :: Rows(:), Cols(:),Diag(:)
2624
    LOGICAL :: Found, Positive
2625
    INTEGER :: pcount
2626
    TYPE(Solver_t), POINTER :: Solver
2627
    TYPE(element_t), POINTER :: Element
2628
    LOGICAL, ALLOCATABLE :: ActiveNodes(:), RowFound(:)
2629

2630
    TYPE(Matrix_t), POINTER :: im
2631

2632
#if defined(ELMER_HAVE_MPIF_HEADER)
2633
    INCLUDE "mpif.h"
2634
#endif
2635

2636
    CALL Info('CRS_FCTLowOrder','Making low order FCT correction to matrix',Level=5)
2637

2638
    N = A % NumberOfRows
2639
    Rows => A % Rows
2640
    Cols => A % Cols
2641
    Diag => A % Diag
2642

2643
    Solver => CurrentModel % Solver
2644
    ALLOCATE(ActiveNodes(n)); activeNodes=.FALSE.
2645

2646
    DO i=1,Solver % NumberofActiveElements
2647
      Element => Solver % Mesh % Elements(Solver % ActiveElements(i))
2648
      IF ( Element % PartIndex /= ParEnv % MyPE ) CYCLE
2649
      ActiveNodes(Solver % Variable % Perm(Element % NodeIndexes)) = .TRUE.
2650
    END do
2651

2652
    IF(.NOT. ASSOCIATED(A % FCT_D) ) THEN
2653
      ALLOCATE( A % FCT_D(SIZE( A % Values) ) ) 
2654
    END IF
2655
    A % FCT_D = 0.0_dp    
2656

2657
    IF(.NOT. ASSOCIATED(A % BulkValues) ) THEN
2658
      ALLOCATE( A % BulkValues(SIZE( A % Values) ) ) 
2659
    END IF
2660
    A % BulkValues = A % Values
2661

2662
    pcount = 0
2663
    Positive = .TRUE.
2664

2665
    N = A % NumberOfRows
2666
    Rows => A % Rows
2667
    Cols => A % Cols
2668
    Diag => A % Diag
2669

2670
    DO i=1,n
2671
      IF ( .NOT. ActiveNodes(i) ) CYCLE
2672

2673
      Aii = A % Values(A % Diag(i))
2674
      IF( Aii > 0.0_dp ) pcount = pcount + 1
2675
      DO k = Rows(i), Rows(i+1)-1
2676
        j = Cols(k)
2677
        IF ( .NOT. ActiveNodes(j) ) CYCLE
2678

2679
        IF( i >= j ) CYCLE
2680

2681
        ! First find entry (j,i)
2682
        Found = .FALSE.
2683
        DO k2 = Rows(j), Rows(j+1)-1
2684
          IF( Cols(k2) == i ) THEN
2685
            Found = .TRUE.
2686
            EXIT
2687
          END IF
2688
        END DO
2689
        IF(.NOT. Found ) THEN
2690
          CALL Fatal('CRS_FCTLowOrder','Entry not found, matrix topology not symmetric!?')
2691
        END IF
2692

2693
        ! Modified so that it also works similarly to A and -A
2694
        Aij = A % Values(k); Aji = A % Values(k2)
2695
        IF (ParEnv % PEs>1) THEN
2696
          Aij = Aij + A % HaloValues(k)
2697
          Aji = Aji + A % HaloValues(k2)
2698
        END IF
2699

2700
        ! Formula (30) in Kuzmin's paper
2701
        ! Positive = Aii > 0.0_dp
2702
        ! In Kuzmin's paper matrix K is -K compared to Elmer convention.
2703
        ! Hence also the condition here is opposite. 
2704
        IF( Positive ) THEN
2705
          Dij = MIN( -Aij, -Aji, 0.0_dp )
2706
        ELSE
2707
          Dij = MAX( -Aij, -Aji, 0.0_dp )
2708
        END IF
2709

2710
        IF(.FALSE.) THEN
2711
          PRINT *,'ij',i,j,Cols(k2),Cols(k)
2712
          PRINT *,'Diag',Cols(Diag(i)),Cols(Diag(j))
2713
          PRINT *,'A',Aij,Aji,Aii,Dij
2714
        END IF
2715

2716
        ! Formula (32) in Kuzmin's paper
2717
        IF( ABS(Dij) > 0._dp ) THEN
2718
          A % FCT_D(k)  = A % FCT_D(k)  + Dij
2719
          A % FCT_D(k2) = A % FCT_D(k2) + Dij
2720
          A % FCT_D(Diag(i)) = A % FCT_D(Diag(i)) - Dij
2721
          A % FCT_D(Diag(j)) = A % FCT_D(Diag(j)) - Dij
2722
        END IF
2723
      END DO
2724
    END DO
2725

2726
    WRITE( Message,'(A,I0,A,I0,A)') 'Positive diagonals ',pcount,' (out of ',n,')'
2727
    CALL Info('CRS_FCTLowOrder',Message,Level=8)
2728

2729
    A % Values = A % Values + A % FCT_D
2730

2731
    ! Just some optional stuff for debugging purposes
2732
    IF (.FALSE.) THEN
2733
      CALL CRS_RowSumInfo( A, A % BulkValues ) 
2734
      CALL CRS_RowSumInfo( A, A % Values ) 
2735
      CALL CRS_RowSumInfo( A, A % FCT_D ) 
2736
    END IF
2737

2738
    ! Create a lumped mass matrix by computing the rowsums of the 
2739
    ! initial mass matrix.
2740
    CALL Info('CRS_FCTLowOrder','Creating lumped mass matrix',Level=10)
2741
    IF(.NOT. ASSOCIATED(A % MassValuesLumped)) THEN         
2742
      ALLOCATE(A % MassValuesLumped(n))
2743
    END IF
2744
    ML => A % MassValuesLumped     
2745
    DO i=1,n
2746
      msum = 0.0_dp
2747
      DO j=Rows(i),Rows(i+1)-1
2748
        msum = msum + A % MassValues(j)
2749
      END DO
2750
      ML(i) = msum 
2751
    END DO
2752

2753
!------------------------------------------------------------------------------
2754
  END SUBROUTINE CRS_FCTLowOrder
2755
!------------------------------------------------------------------------------
2756

2757
 
2758

2759
!------------------------------------------------------------------------------
2760
!> Copies the matrix topology from matrix A to build matrix B. Note that the 
2761
!> topology is really reused, so that if matrix A is destroyed also matrix B 
2762
!> becomes unusable.
2763
!------------------------------------------------------------------------------
2764
  SUBROUTINE CRS_CopyMatrixTopology(A,B)
2765
!------------------------------------------------------------------------------
2766
    TYPE(Matrix_t), INTENT(IN) :: A  !< Initial matrix
2767
    TYPE(Matrix_t) :: B              !< New matrix with same matrix topology
2768
!------------------------------------------------------------------------------
2769
    INTEGER :: kb,n,istat
2770
    N = A % NumberOfRows
2771

2772
    IF ( n == 0 ) THEN
2773
      CALL Fatal('CRS_CopyMatrixTopology','The first matrix is assumed to exist')
2774
    END IF
2775
    
2776
    IF ( A % FORMAT /= MATRIX_CRS ) THEN
2777
      CALL Fatal('CRS_CopyMatrixTopology','The matrix structure should be CRS!')
2778
    END IF
2779
    IF ( B % NumberOfRows /= 0 ) THEN
2780
      CALL Fatal('CRS_CopyMatrixTopology','The other matrix is assumed not to exist')
2781
    END IF
2782

2783
    CALL Info('CRS_CopyMatrixTopology','Reusing matrix topology',Level=9)
2784

2785

2786
    B % NumberOfRows = n
2787
    B % ListMatrix => NULL()
2788
    B % FORMAT = A % FORMAT
2789

2790
    B % Rows => A % Rows
2791
    B % Cols => A % Cols
2792
    IF( ASSOCIATED( A % Diag ) ) THEN
2793
      B % Diag => A % Diag
2794
    END IF
2795

2796
!    IF( ASSOCIATED( A % rhs ) ) THEN
2797
!      ALLOCATE( B % rhs(n), STAT = istat )
2798
!      IF( istat /= 0 ) CALL Fatal('CRS_CopyMatrixTopology','memory allocation error 1')
2799
!      B % rhs = 0.0_dp
2800
!    END IF
2801

2802
    kb = SIZE( A % Values )
2803
    ALLOCATE( B % Values(kb), STAT=istat )
2804
    IF( istat /= 0 ) CALL Fatal('CRS_CopyMatrixTopology','memory allocation error 2')
2805
    B % Values = 0.0_dp
2806

2807
!------------------------------------------------------------------------------
2808
  END SUBROUTINE CRS_CopyMatrixTopology
2809
!------------------------------------------------------------------------------
2810

2811

2812
!------------------------------------------------------------------------------
2813
!> Picks a block from matrix A to build matrix B. It is assumed that the 
2814
!> matrix is split into given number of equally sized blocks.
2815
!------------------------------------------------------------------------------
2816
  SUBROUTINE CRS_BlockMatrixPick(A,B,Blocks,Nrow,Ncol,PickPrec)
2817
!------------------------------------------------------------------------------
2818
    TYPE(Matrix_t), INTENT(IN) :: A   !< Initial matrix
2819
    TYPE(Matrix_t) :: B   !< Submatrix picked from the larger matrix
2820
    INTEGER, INTENT(IN) :: Blocks     !< Number of blocks in the initial matrix
2821
    INTEGER, INTENT(IN) :: Nrow       !< Row to be picked
2822
    INTEGER, INTENT(IN) :: Ncol       !< Column to be picked
2823
    LOGICAL, INTENT(IN), OPTIONAL :: PickPrec
2824
!------------------------------------------------------------------------------    
2825
    INTEGER :: i,j,k,l,kb,n,Nrow0,Ncol0,nsub
2826
    INTEGER :: lsub,isub,istat,modNcol
2827
    LOGICAL :: NewMatrix, Diagonal, DoPrec
2828

2829
    IF(Blocks <= 1) THEN
2830
      CALL Fatal('CRS_BlockMatrixPick','No applicable to just one block!')
2831
      RETURN
2832
    END IF
2833

2834
    CALL Info('CRS_BlockMatrixPick','Picking block ('//I2S(Nrow)//&
2835
        ','//I2S(Ncol)//') from matrix',Level=10)
2836

2837
    DoPrec = .FALSE.
2838
    IF(PRESENT(PickPrec)) DoPrec = PickPrec .AND. ASSOCIATED(A % PrecValues)
2839
    
2840
    
2841
    N = A % NumberOfRows
2842
    Nsub = N / Blocks
2843
    modNcol = MOD( Ncol,Blocks)
2844

2845
    NewMatrix = ( B % NumberOfRows == 0 ) 
2846
    Diagonal = ( Nrow == Ncol ) 
2847

2848
    IF( NewMatrix ) THEN
2849
      CALL Info('CRS_BlockMatrixPick','Allocating new matrix',Level=12)
2850
      B % ListMatrix => NULL()
2851
      B % FORMAT = MATRIX_CRS
2852

2853
      B % NumberOfRows = Nsub    
2854
      kb = 0
2855
      
2856
      DO isub=1,Nsub
2857
        i = Blocks * ( isub - 1 ) + Nrow 
2858
        DO k= A % Rows(i), A % Rows(i+1)-1
2859
          l = A % Cols(k)
2860
          IF( MOD(l,Blocks) == modNcol ) THEN
2861
            kb = kb + 1
2862
          END IF
2863
        END DO
2864
      END DO
2865
      
2866
      IF( kb == 0 ) THEN
2867
        CALL Warn('CRS_BlockMatrixPick','No matrix entries in submatrix')
2868
        RETURN
2869
      END IF
2870

2871
      ALLOCATE(B % Rows(nsub+1),B % Cols(kb), B % Values(kb),STAT=istat )
2872
      IF( istat /= 0 ) CALL Fatal('CRS_BlockMatrixPick','memory allocation error for matrix')
2873

2874
      IF(DoPrec) THEN
2875
        ALLOCATE(B % PrecValues(kb),STAT=istat )
2876
        IF( istat /= 0 ) CALL Fatal('CRS_BlockMatrixPick','memory allocation error for precvalues')
2877
      END IF
2878
    ELSE
2879
      CALL Info('CRS_BlockMatrixPick','Using existing matrix structure',Level=12)
2880
    END IF
2881

2882
    IF( Diagonal ) THEN
2883
      IF( .NOT. ASSOCIATED( B % Diag ) ) THEN
2884
        ALLOCATE( B % Diag(nsub), STAT=istat)
2885
        IF( istat /= 0 ) CALL Fatal('CRS_BlockMatrixPick','memory allocation error for diag')      
2886
      END IF
2887
      IF( .NOT. ASSOCIATED( B % Rhs ) ) THEN
2888
        ALLOCATE( B % rhs(nsub), STAT=istat)
2889
        IF( istat /= 0 ) CALL Fatal('CRS_BlockMatrixPick','memory allocation error rhs')      
2890
      END IF
2891
    END IF
2892

2893

2894
    kb = 1
2895
    DO isub=1,Nsub
2896

2897
      IF( NewMatrix ) B % Rows(isub) = kb
2898
      i = Blocks * ( isub - 1 ) + Nrow 
2899

2900
      DO k = A % Rows(i), A % Rows(i+1)-1
2901

2902
        l = A % Cols(k)
2903
        IF( MOD( l, Blocks ) == modNcol ) THEN
2904
          lsub = ( l - 1) / Blocks + 1
2905

2906
          B % Values(kb) = A % Values(k)
2907
          IF(DoPrec) B % PrecValues(kb) = A % PrecValues(k)
2908
          
2909
          IF( NewMatrix ) THEN
2910
            B % Cols(kb) = lsub
2911
            IF( Diagonal .AND. isub == lsub ) B % Diag(isub) = kb
2912
          END IF
2913
          kb = kb + 1  
2914
        END IF
2915
      END DO
2916
      
2917
      IF( Diagonal ) B % rhs(isub) = A % rhs(i)
2918

2919
    END DO
2920
    IF( NewMatrix ) B % Rows(Nsub+1) = kb
2921

2922
!------------------------------------------------------------------------------
2923
  END SUBROUTINE CRS_BlockMatrixPick
2924
!------------------------------------------------------------------------------
2925

2926

2927
!------------------------------------------------------------------------------
2928
!> Picks a block from matrix A to build matrix B. It is assumed that the 
2929
!> matrix is split by intervals given by the users. For example for AV matrix
2930
!> the user would give the size of V as the input and choose then blocks
2931
!> (1,1), (1,2), (2,1) or (2,2). This logic assumes that nodes are numbered
2932
!> first, followed by other dofs. 
2933
!------------------------------------------------------------------------------
2934
  SUBROUTINE CRS_PartMatrixPick(A,B,Splits,Nrow,Ncol,PreserveColumnIndex)
2935
!------------------------------------------------------------------------------
2936
    TYPE(Matrix_t), INTENT(IN) :: A   !< Initial matrix
2937
    TYPE(Matrix_t), INTENT(OUT) :: B  !< Submatrix picked from the larger matrix
2938
    INTEGER, INTENT(IN) :: Splits(:)  !< Position for the splits
2939
    INTEGER, INTENT(IN) :: Nrow       !< Row to be picked
2940
    INTEGER, INTENT(IN) :: Ncol       !< Column to be picked
2941
    LOGICAL, INTENT(IN) :: PreserveColumnIndex !< Whether to start column index from 1 or preserve it
2942
!------------------------------------------------------------------------------    
2943
    INTEGER :: blocks, i,j,k,l,kb,n,kb0
2944
    INTEGER :: lsub,isub,istat,n1,n2,m1,m2,nsub,msub
2945
    LOGICAL :: NewMatrix, Diagonal
2946
    REAL(KIND=dp) :: PickRatio
2947

2948
    blocks = SIZE( Splits ) + 1
2949
    
2950
    CALL Info('CRS_PartMatrixPick','Picking block ('//I2S(Nrow)//','//I2S(Ncol)//&
2951
        ') part out of ('//I2S(blocks)//','//I2S(blocks)//')',Level=6)
2952

2953
    N = A % NumberOfRows
2954

2955
    IF(blocks <= 1) THEN
2956
      CALL Fatal('CRS_PartMatrixPick','No applicable to just one block!')
2957
    END IF
2958
    IF( Nrow > blocks .OR. Nrow < 1 ) THEN
2959
      CALL Fatal('CRS_PartMatrixPick','Invalid value for Nrow: '//I2S(Nrow))
2960
    END IF
2961
    IF( Ncol > blocks .OR. Ncol < 1) THEN
2962
      CALL Fatal('CRS_PartMatrixPick','Invalid value for Ncol: '//I2S(Nrow))
2963
    END IF
2964

2965
    i = MINVAL( Splits ) 
2966
    IF( i <= 0 ) THEN
2967
      CALL Fatal('CRS_PartMatrixPick','Split must be positive: '//I2S(i))
2968
    END IF
2969
    i = MAXVAL( Splits ) 
2970
    IF( i >= n ) THEN
2971
      CALL Fatal('CRS_PartMatrixPick','Split must be smaller than matrix size: '//I2S(i))
2972
    END IF
2973

2974
    kb0 = A % Rows(n+1) - 1
2975
    CALL Info('CRS_PartMatrixPick','Number of nonzeros in initial matrix: '//I2S(kb0),Level=7)
2976

2977
    IF( Nrow == 1 ) THEN
2978
      n1 = 1 
2979
    ELSE
2980
      n1 = splits(nrow-1) + 1
2981
    END IF
2982
    IF( Nrow == blocks ) THEN
2983
      n2 = n
2984
    ELSE 
2985
      n2 = splits(nrow)
2986
    END IF
2987
    nsub = n2 - n1 + 1
2988
    CALL Info('CRS_PartMatrixPick',&
2989
        'Picking rows from '//I2S(n1)//' to '//I2S(n2),Level=7)
2990
    
2991
    IF( Ncol == 1 ) THEN
2992
      m1 = 1 
2993
    ELSE
2994
      m1 = splits(ncol-1) + 1
2995
    END IF
2996
    IF( Ncol == blocks ) THEN
2997
      m2 = n
2998
    ELSE 
2999
      m2 = splits(ncol)
3000
    END IF   
3001
    msub = m2 - m1 + 1
3002
    CALL Info('CRS_PartMatrixPick',&
3003
        'Picking columns from '//I2S(m1)//' to '//I2S(m2),Level=7)
3004

3005
    CALL Info('CRS_PartMatrixPick',&
3006
        'Sizes of submatrix is '//I2S(nsub)//' x '//I2S(msub),Level=7)
3007

3008

3009
    NewMatrix = ( B % NumberOfRows == 0 ) 
3010
    Diagonal = ( Nrow == Ncol ) 
3011

3012
    IF( NewMatrix ) THEN
3013
      B % ListMatrix => NULL()
3014
      B % FORMAT = MATRIX_CRS
3015

3016
      B % NumberOfRows = Nsub    
3017
      kb = 0
3018
      
3019
      DO i=n1,n2
3020
        DO k= A % Rows(i), A % Rows(i+1)-1
3021
          l = A % Cols(k)
3022
          IF( l >= m1 .AND. l <= m2 ) THEN
3023
            kb = kb + 1
3024
          END IF
3025
        END DO
3026
      END DO
3027
      
3028
      IF( kb == 0 ) THEN
3029
        CALL Warn('CRS_PartMatrixPick','No matrix entries in submatrix')
3030
        RETURN
3031
      END IF
3032

3033
      CALL Info('CRS_PartMatrixPick','Number of nonzeros in submatrix: '//I2S(kb))
3034

3035
      ALLOCATE(B % Rows(nsub+1),B % Cols(kb), B % Values(kb),STAT=istat )
3036
      IF( istat /= 0 ) CALL Fatal('CRS_PartMatrixPick','memory allocation error 1')
3037
    END IF
3038

3039
    IF( Diagonal ) THEN
3040
      IF( .NOT. ASSOCIATED( B % Diag ) ) THEN
3041
        ALLOCATE( B % Diag(nsub), STAT=istat)
3042
        IF( istat /= 0 ) CALL Fatal('CRS_PartkMatrixPick','memory allocation error 2')      
3043
      END IF
3044
      IF( .NOT. ASSOCIATED( B % Rhs ) ) THEN
3045
        ALLOCATE( B % rhs(nsub), STAT=istat)
3046
        IF( istat /= 0 ) CALL Fatal('CRS_PartMatrixPick','memory allocation error 3')      
3047
      END IF
3048
    END IF
3049

3050
    kb = 1
3051
    DO i=n1,n2
3052
      isub = i-n1+1
3053
      IF( NewMatrix ) B % Rows(isub) = kb 
3054
      DO k= A % Rows(i), A % Rows(i+1)-1
3055
        l = A % Cols(k)
3056
        IF( l >= m1 .AND. l <= m2 ) THEN
3057
          B % Values(kb) = A % Values(k)
3058
          IF( NewMatrix ) THEN
3059
            IF( PreserveColumnIndex ) THEN
3060
              lsub = l
3061
            ELSE
3062
              lsub = l-m1+1
3063
            END IF
3064
            B % Cols(kb) = lsub
3065
            IF( Diagonal .AND. isub == lsub ) B % Diag(isub) = kb
3066
          END IF
3067
          kb = kb + 1
3068
        END IF
3069
      END DO      
3070
      IF( Diagonal ) B % rhs(isub) = A % rhs(i)
3071
    END DO
3072
    IF( NewMatrix ) B % Rows(isub+1) = kb
3073
    kb = kb - 1
3074

3075
    PickRatio = 1.0_dp * kb / kb0 
3076
    WRITE( Message,'(A,F8.3,A)') 'Pick matrix ratio is: ',100*PickRatio,' %'
3077
    CALL Info('CRS_PartMatrixPick',Message,Level=6)
3078
    
3079
!------------------------------------------------------------------------------
3080
  END SUBROUTINE CRS_PartMatrixPick
3081
!------------------------------------------------------------------------------
3082

3083

3084
!------------------------------------------------------------------------------
3085
!> Picks a block from matrix A to build matrix B. 
3086
!> This subroutine enables the use of 
3087
!> nontrivial block decompositions. 
3088
!------------------------------------------------------------------------------
3089
  SUBROUTINE CRS_BlockMatrixPick2(A,B,BlockStruct,Nrow,Ncol,PickPrec)
3090
!------------------------------------------------------------------------------
3091
    TYPE(Matrix_t), INTENT(IN) :: A   !< Initial matrix
3092
    TYPE(Matrix_t) :: B   !< Submatrix picked from the larger matrix
3093
    INTEGER, POINTER :: BlockStruct(:)     !< Block decomposition structure of the initial matrix
3094
    INTEGER, INTENT(IN) :: Nrow       !< Row to be picked
3095
    INTEGER, INTENT(IN) :: Ncol       !< Column to be picked
3096
    LOGICAL, INTENT(IN), OPTIONAL :: PickPrec
3097
!------------------------------------------------------------------------------
3098
    INTEGER :: i,j,k,l,kb,n,Nrow0,Ncol0,nsub,Mrow,Mcol,mr,mc,imsub,lmsub
3099
    INTEGER :: lsub,isub,istat,modNcol,Blocks
3100
    LOGICAL :: NewMatrix, Allocated, Diagonal, Hit, DoPrec
3101
    INTEGER, ALLOCATABLE :: Irow(:), Icol(:)
3102
    
3103
    Blocks = SIZE( BlockStruct )
3104

3105
    IF(Blocks <= 1) THEN
3106
      CALL Fatal('CRS_BlockMatrixPick2','Not applicable for just one block!')
3107
      RETURN
3108
    END IF
3109

3110
    DoPrec = .FALSE.
3111
    IF(PRESENT(PickPrec)) DoPrec = PickPrec .AND. ASSOCIATED(A % PrecValues)
3112
    
3113
    N = A % NumberOfRows
3114

3115
    Mrow = 0
3116
    Mcol = 0
3117
    ALLOCATE( Irow(Blocks), Icol(Blocks) ) 
3118
    Irow = 0
3119
    Icol = 0
3120

3121
    DO i=1,Blocks
3122
      IF( BlockStruct(i) == Nrow ) THEN
3123
        Mrow = Mrow + 1
3124
        Irow(Mrow) = i
3125
      END IF
3126
      IF( BlockStruct(i) == Ncol ) THEN
3127
        Mcol = Mcol + 1
3128
        Icol(Mcol) = i
3129
      END IF
3130
    END DO
3131

3132
    IF( Mrow == 0 .OR. Mcol == 0 ) THEN
3133
      CALL Fatal('CRS_BlockMatrixPick2','Nothing to pick!')
3134
    END IF
3135

3136
    Nsub = N / Blocks
3137
    modNcol = MOD( Ncol,Blocks)
3138

3139
    NewMatrix = ( B % NumberOfRows == 0 ) 
3140
    Allocated = .NOT. NewMatrix
3141
    Diagonal = ( Nrow == Ncol ) 
3142

3143
    IF( .NOT. Allocated ) THEN
3144
      !PRINT *,'block rows no:',Mrow,' inds:',Irow(1:Mrow)
3145
      !PRINT *,'block cols no:',Mcol,' inds:',Icol(1:Mcol)
3146
      B % ListMatrix => NULL()
3147
      B % FORMAT = MATRIX_CRS
3148
      B % NumberOfRows = Mrow *  Nsub    
3149
    END IF
3150

3151
100 kb = 1      
3152
    DO isub=1,Nsub
3153

3154
      DO mr=1,Mrow
3155
        imsub = Mrow*(isub-1)+mr
3156
        IF( Allocated .AND. NewMatrix ) B % Rows( imsub ) = kb
3157
        i = Blocks * ( isub - 1 ) + Irow(mr) 
3158
        
3159
        DO k= A % Rows(i), A % Rows(i+1)-1
3160
          l = A % Cols(k)
3161
          Hit = .FALSE.
3162

3163
          DO mc = 1, Mcol
3164
            IF( MOD(l,Blocks) == MOD( Icol(mc), Blocks) ) THEN
3165
              Hit = .TRUE.
3166
              EXIT
3167
            END IF
3168
          END DO
3169

3170
          IF( Hit ) THEN
3171
            IF( Allocated ) THEN
3172
              lmsub = Mcol * ( ( l - 1) / Blocks ) + mc
3173

3174
              B % Values(kb) = A % Values(k)
3175
              IF(DoPrec) B % PrecValues(kb) = A % PrecValues(k)
3176

3177
              IF( NewMatrix ) THEN
3178
                B % Cols(kb) = lmsub
3179
                IF( Diagonal ) THEN
3180
                  IF( imsub == lmsub ) B % Diag(imsub) = kb
3181
                END IF
3182
              END IF            
3183
              IF( Diagonal ) B % rhs(imsub) = A % rhs(i)              
3184
            END IF
3185
            kb = kb + 1
3186
          END IF
3187

3188
        END DO
3189
      END DO
3190
    END DO
3191
    
3192
    IF( .NOT. Allocated ) THEN
3193
      IF( kb == 1 ) THEN
3194
        CALL Warn('CRS_BlockMatrixPick2','No matrix entries in submatrix')
3195
        RETURN
3196
      END IF
3197

3198
      ALLOCATE(B % Rows(Mrow*nsub+1),B % Cols(kb-1), B % Values(kb-1),STAT=istat )
3199
      IF( istat /= 0 ) CALL Fatal('CRS_BlockMatrixPick2','memory allocation error 1')
3200
      
3201
      B % Rows(Mrow*Nsub+1) = kb
3202
      
3203
      IF( Diagonal ) THEN
3204
        ALLOCATE( B % Diag(Mrow*nsub), B % rhs(Mrow*nsub), STAT=istat)
3205
        IF( istat /= 0 ) CALL Fatal('CRS_BlockMatrixPick2','memory allocation error 2')      
3206
      END IF
3207

3208
      IF(DoPrec) THEN
3209
        ALLOCATE(B % PrecValues(kb-1),STAT=istat )
3210
        IF( istat /= 0 ) CALL Fatal('CRS_BlockMatrixPick2','memory allocation error 3')
3211
      END IF
3212
      
3213
      
3214
      IF( A % COMPLEX ) THEN
3215
        IF( MOD( Mrow, 2) == 0 .AND. MOD( Mcol, 2) == 0 ) THEN
3216
          B % COMPLEX = .TRUE.
3217
        END IF
3218
      END IF
3219

3220
      IF( A % Ndeg > 1 ) THEN
3221
        IF( Mrow == Mcol ) THEN
3222
          B % Ndeg = Mrow
3223
        END IF
3224
      END IF
3225
      
3226
      Allocated = .TRUE.
3227
      GOTO 100
3228
    END IF
3229

3230

3231
!------------------------------------------------------------------------------
3232
  END SUBROUTINE CRS_BlockMatrixPick2
3233
!------------------------------------------------------------------------------
3234

3235

3236
!------------------------------------------------------------------------------
3237
!> Copies some preconditioning structures from matrix A to B. 
3238
!> The intent is to allow saving of memory and CPU time for cases
3239
!> where similar preconditioning could be used for many components.
3240
!------------------------------------------------------------------------------
3241
  FUNCTION CRS_CopyMatrixPrec(A,B) RESULT(Status)
3242
!------------------------------------------------------------------------------
3243
    TYPE(Matrix_t), INTENT(IN) :: A  !< Initial matrix
3244
    TYPE(Matrix_t) :: B  !< New matrix with the same preconditioner
3245
    LOGICAL :: Status  !< Returns true if copying was successful or unnecessary
3246
!------------------------------------------------------------------------------
3247
    INTEGER :: n
3248

3249
    Status = .FALSE.
3250

3251
    IF( ASSOCIATED( B % IluValues ) ) THEN
3252
!      CALL Info('CRS_CopyMatrixPrec','ILU preconditioner already exists')
3253
      Status = .TRUE.
3254
    END IF
3255

3256
    IF( ASSOCIATED( B % Ematrix ) ) THEN
3257
!      CALL Info('CRS_CopyMatrixPrec','AMG preconditioner already exists')
3258
      Status = .TRUE.      
3259
    END IF
3260
    
3261
    IF( Status ) RETURN
3262

3263
    IF( SIZE( A % Values ) /= SIZE( B % Values ) ) THEN
3264
      !PRINT *,'sizes',SIZE( A % Values ), SIZE( B % Values )
3265
      CALL Info('CRS_CopyMatrixPrec','Mismatch in size, returning')            
3266
      RETURN
3267
    END IF
3268
    
3269
    IF( ASSOCIATED( A % IluValues ) ) THEN
3270
      CALL Info('CRS_CopyMatrixPrec','Reusing ILU preconditioner topology',Level=9)      
3271
      B % IluRows => A % IluRows
3272
      B % IluCols => A % IluCols
3273
      B % IluDiag => A % IluDiag
3274

3275
      n = SIZE( A % ILUValues ) 
3276
      ALLOCATE( B % IluValues(n) )
3277
      B % IluValues = 0.0_dp
3278
      Status = .TRUE.
3279
      RETURN
3280
    END IF
3281

3282
    RETURN
3283

3284
    ! This should be still worked on....
3285
    !------------------------------------------------
3286
    IF( ASSOCIATED( A % Ematrix ) ) THEN
3287
      CALL Info('CRS_CopyMatrixPrec','Reusing AMG preconditioner topology',Level=9)      
3288
      B % Ematrix => A % Ematrix
3289
    END IF
3290

3291
  END FUNCTION CRS_CopyMatrixPrec
3292

3293
!------------------------------------------------------------------------------
3294
!> Creates CRS matrix of different size but same structure.
3295
!------------------------------------------------------------------------------
3296

3297

3298
  SUBROUTINE CRS_CreateChildMatrix( ParentMat, ParentDofs, ChildMat, Dofs, ColDofs, &
3299
      CreateRhs, NoReuse, Diagonal ) 
3300

3301
    TYPE(Matrix_t) :: ParentMat
3302
    INTEGER :: ParentDofs
3303
    INTEGER :: Dofs
3304
    TYPE(Matrix_t), POINTER :: ChildMat
3305
    INTEGER, OPTIONAL :: ColDofs
3306
    LOGICAL, OPTIONAL :: CreateRhs
3307
    LOGICAL, OPTIONAL :: NoReuse
3308
    LOGICAL, OPTIONAL :: Diagonal
3309

3310
    INTEGER :: i,j,ii,jj,k,l,m,n,nn,Cdofs,Cmult
3311
    LOGICAL :: ReuseMatrix
3312
    LOGICAL :: IsDiagonal
3313
    
3314
    IF( PRESENT( ColDofs ) ) THEN
3315
      CDofs = ColDofs
3316
    ELSE
3317
      CDofs = Dofs
3318
    END IF
3319

3320
    IF( PRESENT( Diagonal ) ) THEN
3321
      IsDiagonal = Diagonal
3322
    ELSE
3323
      IsDiagonal = .FALSE.
3324
    END IF
3325
    
3326
    ReuseMatrix = ( Dofs == ParentDofs .AND. CDofs == ParentDofs )
3327
    IF( PRESENT( NoReuse ) ) THEN
3328
      IF( NoReuse ) ReuseMatrix = .FALSE.         
3329
    END IF
3330

3331

3332
    IF( ReuseMatrix ) THEN
3333
      CALL Info('CRS_CreateChildMatrix','Reusing initial matrix topology',Level=8)    
3334

3335
      ChildMat % Cols => ParentMat % Cols
3336
      ChildMat % Rows => ParentMat % Rows
3337
      ChildMat % Diag => ParentMat % Diag
3338

3339
      ChildMat % NumberOfRows = ParentMat % NumberOfRows
3340

3341
      m = SIZE( ParentMat % Values )
3342
      ALLOCATE( ChildMat % Values(m) )
3343
      ChildMat % Values = 0.0_dp
3344

3345
    ELSE IF( Dofs == ParentDofs .AND. Cdofs == ParentDofs ) THEN
3346
      CALL Info('CRS_CreateChildMatrix','Copying initial matrix topology',Level=8)    
3347

3348
      ALLOCATE( ChildMat % Cols( SIZE(ParentMat % Cols) ) )
3349
      ALLOCATE( ChildMat % Rows( SIZE(ParentMat % Rows) ) )
3350
      ALLOCATE( ChildMat % Diag( SIZE(ParentMat % Diag) ) )
3351

3352
      ChildMat % Cols = ParentMat % Cols
3353
      ChildMat % Rows = ParentMat % Rows
3354
      ChildMat % Diag = ParentMat % Diag
3355

3356
      ChildMat % NumberOfRows = ParentMat % NumberOfRows
3357

3358
      m = SIZE( ParentMat % Values )
3359
      ALLOCATE( ChildMat % Values(m) )
3360
      ChildMat % Values = 0.0_dp
3361
    ELSE IF( IsDiagonal ) THEN
3362

3363
      CALL Info('CRS_CreateChildMatrix','Multiplying initial matrix topology for diagonal system',Level=8)    
3364

3365
      IF( CDofs /= Dofs ) THEN
3366
        CALL Fatal('CRS_CreateChildMatrix','Diagonal matrix must be square matrix!')
3367
      END IF
3368
      
3369
      cmult = Dofs / ParentDofs 
3370
      IF( cmult <= 1 .OR. Dofs /= cmult * ParentDofs ) THEN
3371
        CALL Fatal('CRS_CreateChildMatrix','Diagonal child matrix must be a multiple of parent matrix!')        
3372
      END IF
3373
            
3374
      ALLOCATE( ChildMat % Cols( SIZE(ParentMat % Cols) * cmult ) )
3375
      ALLOCATE( ChildMat % Rows( (SIZE(ParentMat % Rows)-1) * cmult + 1 ) )
3376

3377
      ChildMat % NumberOfRows = ParentMat % NumberOfRows * cmult
3378
      
3379
      ii = 0
3380
      jj = 0
3381
      ChildMat % Rows(1) = 1
3382
      DO i=1, ParentMat % NumberOFRows
3383

3384
        DO k=1,cmult
3385

3386
          ii = ii + 1
3387
          DO j=ParentMat % Rows(i), ParentMat % Rows(i+1)-1
3388
            nn = ParentMat % Cols(j)
3389
            jj = jj + 1            
3390
            ChildMat % Cols(jj) = cmult*(nn-1) + k
3391
          END DO
3392

3393
          ChildMat % Rows(ii+1) = jj+1
3394
        END DO
3395
      END DO
3396
      
3397
      ALLOCATE( ChildMat % Values(jj) )
3398
      ChildMat % Values = 0.0_dp
3399

3400
      IF( Dofs == CDofs ) THEN
3401
        ALLOCATE( ChildMat % Diag( SIZE(ParentMat % Diag) * cmult ) )
3402
        DO i=1,ChildMat % NumberOfRows
3403
          DO j=ChildMat % Rows(i), ChildMat % Rows(i+1)-1
3404
            IF (ChildMat % Cols(j) == i) THEN
3405
              ChildMat % Diag(i) = j
3406
              EXIT
3407
            END IF
3408
          END DO
3409
        END DO
3410
      END IF
3411

3412
    ELSE
3413
      CALL Info('CRS_CreateChildMatrix','Multiplying initial matrix topology',Level=8)    
3414

3415
      ALLOCATE( ChildMat % Cols( SIZE(ParentMat % Cols) * Dofs * CDofs / ParentDofs**2 ) )
3416
      ALLOCATE( ChildMat % Rows( (SIZE(ParentMat % Rows)-1) * Dofs / ParentDofs + 1 ) )
3417

3418
      ChildMat % NumberOfRows = ParentMat % NumberOfRows * Dofs / ParentDofs           
3419

3420
      ii = 0
3421
      jj = 0
3422
      ChildMat % Rows(1) = 1
3423
      DO i=1, ParentMat % NumberOFRows, ParentDOFs
3424
        DO k=1,Dofs
3425
          ii = ii + 1
3426
          DO j=ParentMat % Rows(i), ParentMat % Rows(i+1)-1, ParentDOFs
3427
            nn = (ParentMat % Cols(j)-1) / ParentDofs + 1
3428
            DO l=1,CDofs
3429
              jj = jj + 1
3430
              ChildMat % Cols(jj) = Dofs*(nn-1) + l
3431
            END DO
3432
          END DO
3433
          ChildMat % Rows(ii+1) = jj+1
3434
        END DO
3435
      END DO
3436

3437
      ALLOCATE( ChildMat % Values(jj) )
3438
      ChildMat % Values = 0.0_dp
3439

3440
      IF( Dofs == CDofs ) THEN
3441
        ALLOCATE( ChildMat % Diag( SIZE(ParentMat % Diag) * Dofs / ParentDofs ) )      
3442
        DO i=1,ChildMat % NumberOfRows
3443
          DO j=ChildMat % Rows(i), ChildMat % Rows(i+1)-1
3444
            IF (ChildMat % Cols(j) == i) THEN
3445
              ChildMat % Diag(i) = j
3446
              EXIT
3447
            END IF
3448
          END DO
3449
        END DO
3450
      END IF
3451
    END IF
3452

3453
    IF( PRESENT( CreateRhs ) ) THEN
3454
      IF( CreateRhs ) THEN
3455
        ALLOCATE( ChildMat % rhs(ChildMat % NumberOfRows ) )
3456
        ChildMat % rhs = 0.0_dp
3457
      END IF
3458
    END IF
3459

3460
    CALL Info('CRS_CreateChildMatrix','Created matrix with rows: '&
3461
        //I2S( ChildMat % NumberOfRows),Level=10 )
3462

3463

3464
  END SUBROUTINE CRS_CreateChildMatrix
3465
  
3466
  
3467

3468
!------------------------------------------------------------------------------
3469
!>    Builds an incomplete (ILU(n)) factorization for a iterative solver
3470
!>    preconditioner. Real matrix version.
3471
!------------------------------------------------------------------------------
3472
  FUNCTION CRS_IncompleteLU(A,ILUn,Params) RESULT(Status)
3473
!------------------------------------------------------------------------------
3474
    TYPE(Matrix_t) :: A          !< Structure holding input matrix, will also hold the factorization on exit.
3475
    INTEGER, INTENT(IN) :: ILUn  !< Order of fills allowed 0-9
3476
    TYPE(ValueList_t), POINTER, INTENT(in) :: Params !< 
3477
    LOGICAL :: Status            !< Whether or not the factorization succeeded.
3478
!------------------------------------------------------------------------------
3479
    LOGICAL :: Warned, Retry, Found
3480
    INTEGER :: i,j,k,l,m,n,istat
3481
    INTEGER, POINTER :: Cols(:),Rows(:),Diag(:)
3482
    REAL(KIND=dp), POINTER :: ILUValues(:), Values(:)
3483
    REAL(KIND=dp) :: st, tx, scl, cFactor
3484
    LOGICAL, ALLOCATABLE :: C(:), D(:)
3485
    REAL(KIND=dp), ALLOCATABLE ::  S(:), T(:)
3486
    INTEGER, POINTER :: ILUCols(:),ILURows(:),ILUDiag(:)
3487
    TYPE(Matrix_t), POINTER :: A1
3488
!------------------------------------------------------------------------------
3489
    WRITE(Message,'(a,i1,a)')  &
3490
         'ILU(',ILUn,') (Real), Performing Factorization:'
3491
    CALL Info( 'CRS_IncompleteLU', Message, Level = 6 )
3492
    st = CPUTime()
3493

3494
    N = A % NumberOfRows
3495

3496
    IF(N == 0) THEN
3497
       A % ILURows => A % Rows
3498
       A % ILUCols => A % Cols
3499
       A % ILUDiag => A % Diag
3500

3501
       Status = .TRUE.
3502
       RETURN
3503
    END IF
3504

3505
    Diag   => A % Diag
3506
    Rows   => A % Rows
3507
    Cols   => A % Cols
3508
    IF(ASSOCIATED(A % PrecValues)) THEN
3509
      Values => A % PrecValues
3510
    ELSE
3511
      Values => A % Values
3512
    END IF
3513
    
3514
    IF ( .NOT. ASSOCIATED(A % ILUValues) ) THEN
3515

3516
       IF ( ILUn == 0 ) THEN
3517
          A % ILURows => A % Rows
3518
          A % ILUCols => A % Cols
3519
          A % ILUDiag => A % Diag
3520
       ELSE
3521
          CALL InitializeILU1( A,n )
3522

3523
          IF ( ILUn > 1 ) THEN
3524
             ALLOCATE( A1 )
3525

3526
             DO i=1,ILUn-1
3527
                CALL Info('CRS_IncompleLU','Recursive round: '//I2S(i),Level=7)
3528

3529
                A1 % Cols => A % ILUCols
3530
                A1 % Rows => A % ILURows
3531
                A1 % Diag => A % ILUDiag
3532

3533
                CALL InitializeILU1( A1,n )
3534

3535
                A % ILUCols => A1 % ILUCols
3536
                A % ILURows => A1 % ILURows
3537
                A % ILUDiag => A1 % ILUDiag
3538

3539
                DEALLOCATE( A1 % Cols, A1 % Rows, A1 % Diag )
3540
             END DO
3541

3542
             DEALLOCATE(A1)
3543
          END IF
3544
       END IF
3545

3546
       m = A % ILURows(N+1)-1
3547
       ALLOCATE( A % ILUValues(A % ILURows(N+1)-1), STAT=istat )
3548

3549
       IF ( istat /= 0 ) THEN
3550
         CALL Fatal( 'CRS_IncompleteLU', 'Memory allocation error.' )
3551
       ELSE
3552
         CALL Info('CRS_IncompleteLU','Allocated LU matrix of size: '//i2s(m),Level=10 ) 
3553
       END IF
3554
    END IF
3555

3556
    ILURows   => A % ILURows
3557
    ILUCols   => A % ILUCols
3558
    ILUDiag   => A % ILUDiag
3559
    ILUValues => A % ILUValues
3560
!
3561
!   Allocate space for storing one full row:
3562
!   ----------------------------------------
3563
    ALLOCATE( C(n), S(n) )
3564
    C = .FALSE.
3565
    S =  0.0d0
3566

3567
    IF ( A % Cholesky ) THEN
3568
      CALL Info('CRS_IncompleteLU','Performing incomplete Cholesky',Level=12)
3569

3570
      ALLOCATE( T(n) )
3571

3572
      Retry = ListGetLogical( Params, 'Linear System Cholesky Retry', Found )
3573
      IF( Retry ) THEN
3574
        cFactor = ListGetCReal( Params, 'Linear System Cholesky Factor', Found )
3575
        cFactor = cFactor + 1._dp
3576
      END IF
3577

3578
      Scl = 1._dp
3579

3580
1     CONTINUE
3581

3582
     T =  0._dp
3583
     !
3584
     ! The factorization row by row:
3585
     ! -----------------------------
3586
     Warned = .FALSE.
3587
     DO i=1,n
3588

3589
       ! Convert current row to full form for speed,
3590
       ! only flagging the nonzero entries:
3591
       ! -------------------------------------------
3592
       DO k=Rows(i), Diag(i)
3593
         j = Cols(k)
3594
         T(j) = Values(k)
3595
       END DO
3596

3597
       DO k=ILURows(i), ILUDiag(i)
3598
         j = ILUCols(k)
3599
         C(j) = .TRUE.
3600
         S(j) = ILUValues(k)
3601
       END DO
3602

3603
       ! This is the factorization part for the current row:
3604
       ! ---------------------------------------------------
3605
       S(i) = Scl * T(i)
3606
       DO m=ILURows(i),ILUDiag(i)-1
3607
         j = ILUCols(m)
3608
         S(j) = T(j)
3609
         DO l = ILURows(j),ILUDiag(j)-1
3610
           k = ILUCols(l)
3611
           S(j) = S(j) - S(k) * ILUValues(l)
3612
         END DO
3613
         S(j) = S(j) * ILUValues(ILUDiag(j))
3614
         S(i) = S(i) - S(j)**2
3615
       END DO
3616

3617
       IF ( S(i) <= AEPS ) THEN
3618
         S(i) = 1._dp
3619
         IF ( .NOT. Warned )  THEN
3620
           CALL Warn( 'Cholesky factorization:', &
3621
               'Negative diagonal: not pos.def. or badly conditioned matrix' )
3622
           Warned = .TRUE.
3623
         END IF
3624
         IF ( Retry ) THEN
3625
           Scl = Scl * cFactor
3626
           WRITE(Message, *) Scl
3627
           CALL Warn( 'Cholesky factorization:', &
3628
                  'Retry using diagonal scaling:'//TRIM(Message) )
3629
           GOTO 1
3630
         END IF
3631
       ELSE
3632
         S(i) = 1._dp / SQRT(S(i))
3633
       END IF
3634

3635
       ! Convert the row back to  CRS format:
3636
       ! ------------------------------------
3637
       DO k=Rows(i), Diag(i)
3638
         j = Cols(k) 
3639
         T(j) = 0._dp
3640
       END DO
3641

3642
       DO k=ILURows(i), ILUDiag(i)
3643
         j = ILUCols(k)
3644
         ILUValues(k) = S(j)
3645
         S(j) =  0._dp
3646
         C(j) = .FALSE.
3647
       END DO
3648
     END DO
3649

3650
    ELSE
3651
      CALL Info('CRS_IncompleteLU','Performing incomplete LU',Level=12)
3652

3653
     ! The factorization row by row:
3654
     ! -----------------------------
3655
     DO i=1,N
3656

3657
       ! Convert current row to full form for speed,
3658
       ! only flagging the nonzero entries:
3659
       ! -------------------------------------------
3660
       DO k=Rows(i), Rows(i+1)-1
3661
          S(Cols(k)) = Values(k)
3662
       END DO
3663

3664
       DO k = ILURows(i), ILURows(i+1)-1
3665
          C(ILUCols(k)) = .TRUE.
3666
       END DO
3667
!
3668
!      This is the factorization part for the current row:
3669
!      ---------------------------------------------------
3670
       DO m=ILURows(i),ILUDiag(i)-1
3671
         k = ILUCols(m)
3672
         IF ( S(k) == 0._dp ) CYCLE
3673

3674
         IF ( ABS(ILUValues(ILUDiag(k))) > AEPS ) &
3675
           S(k) = S(k) / ILUValues(ILUDiag(k)) 
3676

3677
         DO l = ILUDiag(k)+1, ILURows(k+1)-1
3678
           j = ILUCols(l)
3679
           IF ( C(j) ) THEN
3680
             S(j) = S(j) - S(k) * ILUValues(l)
3681
           END IF
3682
         END DO
3683
       END DO
3684

3685
       
3686
!
3687
!      Convert the row back to  CRS format:
3688
!      ------------------------------------
3689
       DO k=ILURows(i), ILURows(i+1)-1
3690
         IF ( C(ILUCols(k)) ) THEN
3691
           ILUValues(k)  = S(ILUCols(k))
3692
           S(ILUCols(k)) =  0.0d0
3693
           C(ILUCols(k)) = .FALSE.
3694
         END IF
3695
       END DO
3696
     END DO
3697

3698
     ! Prescale the diagonal for the LU solve:
3699
     ! ---------------------------------------
3700
     DO i=1,N
3701
       IF ( ABS(ILUValues(ILUDiag(i))) < AEPS ) THEN
3702
         ILUValues(ILUDiag(i)) = 1.0d0
3703
       ELSE
3704
         ILUValues(ILUDiag(i)) = 1.0d0 / ILUValues(ILUDiag(i))
3705
       END IF
3706
     END DO
3707
    END IF
3708

3709

3710
!------------------------------------------------------------------------------
3711
    IF( InfoActive(20) ) THEN
3712
      PRINT *,'ILU range:',MINVAL(ILUValues),MAXVAL(ILUValues),&
3713
          SUM(ILUValues)/SIZE(ILUValues),SUM(ABS(ILUValues))/SIZE(ILUValues)
3714
    END IF
3715

3716
    WRITE(Message,'(a,i1,a,i9)') 'ILU(', ILUn, &
3717
        ') (Real), NOF nonzeros: ',ILURows(n+1)
3718
    CALL Info( 'CRS_IncompleteLU', Message, Level=6 )
3719

3720
    WRITE(Message,'(a,i1,a,i9)') 'ILU(', ILUn, &
3721
        ') (Real), filling (%) : ',   &
3722
         FLOOR(ILURows(n+1)*(100.0d0/Rows(n+1)))
3723
    CALL Info( 'CRS_IncompleteLU', Message, Level=6 )
3724

3725
    WRITE(Message,'(A,I1,A,F8.2)') 'ILU(',ILUn, &
3726
        ') (Real), Factorization time (s): ', CPUTime()-st
3727
    CALL Info( 'CRS_IncompleteLU', Message, Level=6 )
3728

3729
    Status = .TRUE.
3730
!------------------------------------------------------------------------------
3731

3732
  CONTAINS
3733

3734
!------------------------------------------------------------------------------
3735
    SUBROUTINE InitializeILU1( A, n )
3736
!------------------------------------------------------------------------------
3737
      TYPE(Matrix_t) :: A
3738
      INTEGER :: n
3739
      INTEGER :: i,j,k,l,istat,RowMin,RowMax,Nonzeros
3740
      INTEGER :: C(n)
3741
      INTEGER, POINTER :: Cols(:),Rows(:),Diag(:), &
3742
           ILUCols(:),ILURows(:),ILUDiag(:)
3743
!------------------------------------------------------------------------------
3744

3745
      Diag => A % Diag
3746
      Rows => A % Rows
3747
      Cols => A % Cols
3748

3749
      ALLOCATE( A % ILURows(N+1),A % ILUDiag(N),STAT=istat )
3750
      IF ( istat /= 0 ) THEN
3751
         CALL Fatal( 'CRS_IncompleteLU', 'Memory allocation error.' )
3752
      END IF
3753

3754
      ILURows => A % ILURows
3755
      ILUDiag => A % ILUDiag
3756
!
3757
!     Count fills, row by row:
3758
!     ------------------------
3759
      NonZeros = Rows(N+1) - 1
3760

3761
      C = 0
3762
      DO i=1,n
3763
         DO k=Rows(i), Rows(i+1)-1
3764
            C(Cols(k)) = 1
3765
         END DO
3766

3767
         DO k = Cols(Rows(i)), i-1
3768
            IF ( C(k) /= 0 ) THEN
3769
               DO l=Diag(k)+1, Rows(k+1)-1
3770
                  j = Cols(l)
3771
                  IF ( C(j) == 0 ) THEN
3772
                    Nonzeros = Nonzeros + 1
3773

3774
                    IF( Nonzeros == HUGE( NonZeros ) ) THEN
3775
                      CALL Error('CRS_IncompleteLU','Number of nonzeros larger than HUGE(Integer)')
3776
                      CALL Fatal('CRS_IncompleteLU','Try some cheaper preconditioner!')                     
3777
                    END IF
3778

3779
                  END IF
3780
               END DO
3781
            END IF
3782
         END DO
3783

3784
         DO k = Rows(i), Rows(i+1)-1
3785
            C(Cols(k)) = 0
3786
         END DO
3787
      END DO
3788

3789
      CALL Info('CRS_IncompleteLU','Number of nonzeros: '//I2S(NonZeros),Level=12)
3790

3791
!------------------------------------------------------------------------------
3792

3793
      ALLOCATE( A % ILUCols(Nonzeros),STAT=istat )
3794
      IF ( istat /= 0 ) THEN
3795
         CALL Fatal( 'CRS_IncompleteLU', 'Memory allocation error.' )
3796
      END IF
3797
      ILUCols => A % ILUCols
3798

3799
!------------------------------------------------------------------------------
3800

3801
!
3802
!     Update row nonzero structures: 
3803
!     ------------------------------
3804
      C = 0
3805
      ILURows(1) = 1
3806
      DO i=1,n
3807
         DO k=Rows(i), Rows(i+1)-1
3808
            C(Cols(k)) = 1
3809
         END DO
3810

3811
         RowMin = Cols(Rows(i))
3812
         RowMax = Cols(Rows(i+1)-1)
3813

3814
         DO k=RowMin, i-1
3815
            IF ( C(k) == 1 ) THEN
3816
               DO l=Diag(k)+1,Rows(k+1)-1
3817
                  j = Cols(l)
3818
                  IF ( C(j) == 0 ) THEN
3819
                     C(j) = 2
3820
                     RowMax = MAX( RowMax, j )
3821
                  END IF
3822
               END DO
3823
            END IF
3824
         END DO
3825

3826
         j = ILURows(i) - 1
3827
         DO k = RowMin, RowMax 
3828
            IF ( C(k) > 0 ) THEN
3829
               j = j + 1
3830
               C(k) = 0
3831
               ILUCols(j) = k
3832
               IF ( k == i ) ILUDiag(i) = j
3833
            END IF
3834
         END DO
3835
         ILURows(i+1) = j + 1
3836
      END DO
3837

3838
      CALL Info('CRS_IncompleteLU','Updated nonzero elements',Level=20)
3839

3840
!------------------------------------------------------------------------------
3841
    END SUBROUTINE InitializeILU1
3842
!------------------------------------------------------------------------------
3843
!------------------------------------------------------------------------------
3844
  END FUNCTION CRS_IncompleteLU
3845
!------------------------------------------------------------------------------
3846

3847
!------------------------------------------------------------------------------
3848
!>    Buids an incomplete (ILU(n)) factorization for an iterative solver
3849
!>    preconditioner. Complex matrix version.
3850
!------------------------------------------------------------------------------
3851
  FUNCTION CRS_ComplexIncompleteLU(A,ILUn) RESULT(Status)
3852
!------------------------------------------------------------------------------
3853
    TYPE(Matrix_t) :: A  !< Structure holding input matrix, will also hold the factorization on exit.
3854
    INTEGER, INTENT(IN) :: ILUn   !< Order of fills allowed 0-9
3855
    LOGICAL :: Status  !< Whether or not the factorization succeeded.
3856
!------------------------------------------------------------------------------
3857
    INTEGER :: i,j,k,l,m,n,istat
3858
    INTEGER, POINTER :: Cols(:),Rows(:),Diag(:)
3859
    REAL(KIND=dp), POINTER ::  Values(:)
3860
    COMPLEX(KIND=dp), POINTER :: ILUValues(:)
3861
    INTEGER, POINTER :: ILUCols(:),ILURows(:),ILUDiag(:)
3862
    TYPE(Matrix_t), POINTER :: A1
3863
    REAL(KIND=dp) :: st
3864
    LOGICAL, ALLOCATABLE :: C(:)
3865
    COMPLEX(KIND=dp), ALLOCATABLE :: S(:), T(:)
3866
!------------------------------------------------------------------------------
3867

3868
    WRITE(Message,'(a,i1,a)') 'ILU(',ILUn,') (Complex), Performing Factorization:'
3869
    CALL Info( 'CRS_ComplexIncompleteLU', Message, Level=6 )
3870
    st = CPUTime()
3871

3872
    N = A % NumberOfRows
3873
    Diag   => A % Diag
3874
    Rows   => A % Rows
3875
    Cols   => A % Cols
3876
    IF(ASSOCIATED(A % PrecValues)) THEN
3877
      CALL Info( 'CRS_ComplexIncompleteLU', 'Factorizing PrecValues', Level=20 )
3878
      Values => A % PrecValues
3879
    ELSE
3880
      CALL Info( 'CRS_ComplexIncompleteLU', 'Factorizing the primary matrix', Level=20 )
3881
      Values => A % Values
3882
    END IF
3883

3884
    IF ( .NOT.ASSOCIATED(A % CILUValues) ) THEN
3885

3886
       ALLOCATE( A1 )
3887
       A1 % NumberOfRows = N / 2
3888

3889
       ALLOCATE( A1 % Rows(n/2+1) )
3890
       ALLOCATE( A1 % Diag(n/2) )
3891
       ALLOCATE( A1 % Cols(SIZE(A % Cols) / 4) )
3892

3893
       A1 % Rows(1) = 1
3894
       k = 0
3895
       DO i=1,n,2
3896
          DO j=A % Rows(i),A % Rows(i+1)-1,2
3897
             k = k + 1
3898
             A1 % Cols(k) = (A % Cols(j)+1) / 2
3899
             IF ( A % Cols(j) == i ) A1 % Diag((i+1)/2) = k
3900
          END DO
3901
          A1 % Rows((i+1)/2+1) = k+1
3902
       END DO
3903

3904
       IF ( ILUn == 0 ) THEN
3905
          A % ILUCols => A1 % Cols
3906
          A % ILURows => A1 % Rows
3907
          A % ILUDiag => A1 % Diag
3908
       ELSE
3909
          CALL InitializeComplexILU1( A1, n/2 )
3910

3911
          A % ILUCols => A1 % ILUCols
3912
          A % ILURows => A1 % ILURows
3913
          A % ILUDiag => A1 % ILUDiag
3914

3915
          DEALLOCATE( A1 % Cols,A1 % Rows,A1 % Diag )
3916
       END IF
3917

3918
       DEALLOCATE( A1 )
3919

3920
       IF ( ILUn > 1 ) THEN
3921
          ALLOCATE( A1 )
3922
          A1 % NumberOfRows = N / 2
3923

3924
          DO i=1,ILUn-1
3925
             A1 % Cols => A % ILUCols
3926
             A1 % Rows => A % ILURows
3927
             A1 % Diag => A % ILUDiag
3928

3929
             CALL InitializeComplexILU1( A1, n/2 )
3930

3931
             A % ILUCols => A1 % ILUCols
3932
             A % ILURows => A1 % ILURows
3933
             A % ILUDiag => A1 % ILUDiag
3934

3935
             DEALLOCATE( A1 % Cols,A1 % Rows,A1 % Diag )
3936
          END DO
3937

3938
          DEALLOCATE(A1)
3939
       END IF
3940

3941
       ALLOCATE( A % CILUValues(A % ILURows(N/2+1)),STAT=istat )
3942

3943
       IF ( istat /= 0 ) THEN
3944
          CALL Fatal( 'CRS_ComplexIncompleteLU', 'Memory allocation error.' )
3945
       END IF
3946
    END IF
3947

3948
    ILURows   => A % ILURows
3949
    ILUCols   => A % ILUCols
3950
    ILUDiag   => A % ILUDiag
3951
    ILUValues => A % CILUValues
3952

3953
!
3954
!   Allocate space for storing one full row:
3955
!   ----------------------------------------
3956
    ALLOCATE( C(n/2), S(n/2) )
3957
    C = .FALSE.
3958
    S =  0.0d0
3959

3960
    IF ( A % Cholesky ) THEN
3961
      ALLOCATE( T(n/2) )
3962
      T =  0.0d0
3963
     !
3964
     ! The factorization row by row:
3965
     ! -----------------------------
3966
     DO i=1,N/2
3967

3968
       ! Convert current row to full form for speed,
3969
       ! only flagging the nonzero entries:
3970
       ! -------------------------------------------
3971
       DO k = Rows(2*i-1), Rows(2*i)-1,2
3972
          T((Cols(k)+1)/2) = CMPLX( Values(k), -Values(k+1), KIND=dp )
3973
       END DO
3974

3975
       DO j=ILURows(i), ILUDiag(i)
3976
          C(ILUCols(j)) = .TRUE.
3977
          S(ILUCols(j)) = ILUValues(j)
3978
       END DO
3979

3980
       ! This is the factorization part for the current row:
3981
       ! ---------------------------------------------------
3982
       S(i) = T(i)
3983
       DO m=ILURows(i),ILUDiag(i)-1
3984
         j = ILUCols(m)
3985
         S(j) = T(j)
3986
         DO l = ILURows(j),ILUDiag(j)-1
3987
           k = ILUCols(l)
3988
           S(j) = S(j) - S(k) * CONJG(ILUValues(l))
3989
         END DO
3990
         S(j) = S(j) * ILUValues(ILUDiag(j))
3991
         S(i) = S(i) - S(j)*CONJG(S(j))
3992
       END DO
3993

3994
       S(i) = 1._dp / SQRT(S(i))
3995

3996
       ! Convert the row back to  CRS format:
3997
       ! ------------------------------------
3998
       DO k = Rows(2*i-1), Rows(2*i)-1,2
3999
         T((Cols(k)+1)/2) =  0._dp
4000
       END DO
4001

4002
       DO k=ILURows(i), ILUDiag(i)
4003
         ILUValues(k)  = S(ILUCols(k))
4004
         S(ILUCols(k)) =  0._dp
4005
         C(ILUCols(k)) = .FALSE.
4006
       END DO
4007
     END DO
4008

4009
    ELSE
4010

4011
     ! The factorization row by row:
4012
     ! -----------------------------
4013
     DO i=1,N/2
4014

4015
       ! Convert the current row to full form for speed,
4016
       ! only flagging the nonzero entries:
4017
       ! -----------------------------------------------
4018
       DO k = ILURows(i), ILURows(i+1)-1
4019
         C(ILUCols(k)) = .TRUE.
4020
       END DO
4021

4022
       DO k = Rows(2*i-1), Rows(2*i)-1,2
4023
         S((Cols(k)+1)/2) = CMPLX( Values(k), -Values(k+1), KIND=dp )
4024
       END DO
4025

4026
       ! This is the factorization part for the current row:
4027
       ! ---------------------------------------------------
4028
       DO m=ILURows(i),ILUDiag(i)-1
4029
         k = ILUCols(m)
4030
         IF ( S(k) == 0._dp ) CYCLE
4031

4032
         IF ( ABS(ILUValues(ILUDiag(k))) > AEPS ) &
4033
           S(k) = S(k) / ILUValues(ILUDiag(k)) 
4034

4035
         DO l = ILUDiag(k)+1, ILURows(k+1)-1
4036
           j = ILUCols(l)
4037
           IF ( C(j) ) THEN
4038
             S(j) = S(j) - S(k) * ILUValues(l)
4039
           END IF
4040
         END DO
4041
       END DO
4042

4043
       ! Convert the row back to  CRS format:
4044
       ! ------------------------------------
4045
       DO k=ILURows(i), ILURows(i+1)-1
4046
         ILUValues(k)  = S(ILUCols(k))
4047
         S(ILUCols(k)) =  0._dp
4048
         C(ILUCols(k)) = .FALSE.
4049
       END DO
4050
     END DO
4051

4052
     ! Prescale the diagonal for the LU solve:
4053
     ! ---------------------------------------
4054
     DO i=1,n/2
4055
       IF ( ABS(ILUValues(ILUDiag(i))) < AEPS ) THEN
4056
         ILUValues(ILUDiag(i)) = 1._dp
4057
       ELSE
4058
         ILUValues(ILUDiag(i)) = 1._dp / ILUValues(ILUDiag(i))
4059
       END IF
4060
     END DO
4061
    END IF
4062

4063
!------------------------------------------------------------------------------
4064

4065
    WRITE(Message,'(a,i1,a,i9)') 'ILU(', ILUn, &
4066
        ') (Complex), NOF nonzeros: ',ILURows(n/2+1)
4067
    CALL Info( 'CRS_ComplexIncompleteLU', Message, Level=6 )
4068

4069
    WRITE(Message,'(a,i1,a,i9)') 'ILU(', ILUn, &
4070
        ') (Complex), filling (%) : ',   &
4071
         FLOOR(ILURows(n/2+1)*(400.0d0/Rows(n+1)))
4072
    CALL Info( 'CRS_ComplexIncompleteLU', Message, Level=6 )
4073

4074
    WRITE(Message,'(A,I1,A,F8.2)') 'ILU(',ILUn, &
4075
        ') (Complex), Factorization ready at (s): ', CPUTime()-st
4076
    CALL Info( 'CRS_ComplexIncompleteLU', Message, Level=6 )
4077

4078
    Status = .TRUE.
4079
!------------------------------------------------------------------------------
4080

4081
  CONTAINS
4082

4083
!------------------------------------------------------------------------------
4084
    SUBROUTINE InitializeComplexILU1( A, n )
4085
!------------------------------------------------------------------------------
4086
      TYPE(Matrix_t) :: A
4087
      INTEGER :: n
4088

4089
      INTEGER :: i,j,k,l,istat,RowMin,RowMax,Nonzeros
4090

4091
      INTEGER :: C(n)
4092
      INTEGER, POINTER :: Cols(:),Rows(:),Diag(:), &
4093
           ILUCols(:),ILURows(:),ILUDiag(:)
4094
!------------------------------------------------------------------------------
4095

4096
      Diag => A % Diag
4097
      Rows => A % Rows
4098
      Cols => A % Cols
4099

4100
      ALLOCATE( A % ILURows(N+1),A % ILUDiag(N),STAT=istat )
4101
      IF ( istat /= 0 ) THEN
4102
         CALL Fatal( 'CRS_ComplexIncompleteLU', 'Memory allocation error.' )
4103
      END IF
4104

4105
      ILURows => A % ILURows
4106
      ILUDiag => A % ILUDiag
4107

4108
!
4109
!     Count fills, row by row:
4110
!     ------------------------
4111
      NonZeros = Rows(N+1) - 1
4112
      C = 0
4113
      DO i=1,n
4114
         DO k=Rows(i), Rows(i+1)-1
4115
            C(Cols(k)) = 1
4116
         END DO
4117

4118
         DO k = Cols(Rows(i)), i-1
4119
            IF ( C(k) /= 0 ) THEN
4120
               DO l=Diag(k)+1, Rows(k+1)-1
4121
                  j = Cols(l)
4122
                  IF ( C(j) == 0 ) Nonzeros = Nonzeros + 1
4123
               END DO
4124
            END IF
4125
         END DO
4126

4127
         DO k = Rows(i), Rows(i+1)-1
4128
            C(Cols(k)) = 0
4129
         END DO
4130
      END DO
4131

4132
!------------------------------------------------------------------------------
4133

4134
      ALLOCATE( A % ILUCols(Nonzeros),STAT=istat )
4135
      IF ( istat /= 0 ) THEN
4136
         CALL Fatal( 'CRS_ComplexIncompleteLU', 'Memory allocation error.' )
4137
      END IF
4138
      ILUCols => A % ILUCols
4139

4140
!------------------------------------------------------------------------------
4141

4142
!
4143
!     Update row nonzero structures: 
4144
!     ------------------------------
4145
      C = 0
4146
      ILURows(1) = 1
4147
      DO i=1,n
4148
         DO k=Rows(i), Rows(i+1)-1
4149
            C(Cols(k)) = 1
4150
         END DO
4151

4152
         RowMin = Cols(Rows(i))
4153
         RowMax = Cols(Rows(i+1)-1)
4154

4155
         DO k=RowMin, i-1
4156
            IF ( C(k) == 1 ) THEN
4157
               DO l=Diag(k)+1,Rows(k+1)-1
4158
                  j = Cols(l)
4159
                  IF ( C(j) == 0 ) THEN
4160
                     C(j) = 2
4161
                     RowMax = MAX( RowMax, j )
4162
                  END IF
4163
               END DO
4164
            END IF
4165
         END DO
4166

4167
         j = ILURows(i) - 1
4168
         DO k = RowMin, RowMax 
4169
            IF ( C(k) > 0 ) THEN
4170
               j = j + 1
4171
               C(k) = 0
4172
               ILUCols(j) = k
4173
               IF ( k == i ) ILUDiag(i) = j
4174
            END IF
4175
         END DO
4176
         ILURows(i+1) = j + 1
4177
      END DO
4178
!------------------------------------------------------------------------------
4179
    END SUBROUTINE InitializeComplexILU1
4180
!------------------------------------------------------------------------------
4181
!------------------------------------------------------------------------------
4182
  END FUNCTION CRS_ComplexIncompleteLU
4183
!------------------------------------------------------------------------------
4184

4185

4186
!------------------------------------------------------------------------------
4187
!>    Buids an incomplete (ILUT) factorization for an iterative solver
4188
!>    preconditioner. Real matrix version.
4189
!------------------------------------------------------------------------------
4190
  FUNCTION CRS_ILUT(A,TOL) RESULT(Status)
4191
!------------------------------------------------------------------------------
4192
    TYPE(Matrix_t) :: A  !< Structure holding input matrix, will also hold the factorization on exit.
4193
    REAL(KIND=dp), INTENT(IN) :: TOL  !< Drop toleranece: if ILUT(i,j) <= NORM(A(i,:))*TOL the value is dropped.
4194
    LOGICAL :: Status    !< Whether or not the factorization succeeded.
4195
!------------------------------------------------------------------------------
4196
    INTEGER :: n
4197
    REAL(KIND=dp) :: t
4198
!------------------------------------------------------------------------------
4199

4200
    CALL Info( 'CRS_ILUT', 'Performing factorization:', Level=6 )
4201
    t = CPUTime()
4202

4203
    n = A % NumberOfRows
4204

4205
    IF ( ASSOCIATED( A % ILUValues ) ) THEN
4206
       DEALLOCATE( A % ILURows, A % ILUDiag, A % ILUCols, A % ILUValues )
4207
    END IF
4208
!
4209
!   ... and then to the point:
4210
!   --------------------------
4211
    CALL ComputeILUT( A, n, TOL )
4212
! 
4213
    WRITE( Message, * ) 'ILU(T) (Real), NOF nonzeros: ',A % ILURows(N+1)
4214
    CALL Info( 'CRS_ILUT', Message, Level=6 )
4215
    WRITE( Message, * ) 'ILU(T) (Real), filling (%): ', &
4216
         FLOOR(A % ILURows(N+1)*(100.0d0/A % Rows(N+1)))
4217
    CALL Info( 'CRS_ILUT', Message, Level=6 )
4218
    WRITE(Message,'(A,F8.2)') 'ILU(T) (Real), Factorization ready at (s): ', CPUTime()-t
4219
    CALL Info( 'CRS_ILUT', Message, Level=6 )
4220

4221
    Status = .TRUE.
4222
!------------------------------------------------------------------------------
4223

4224
  CONTAINS
4225

4226
!------------------------------------------------------------------------------
4227
    SUBROUTINE ComputeILUT( A,n,TOL )
4228
!------------------------------------------------------------------------------
4229
      REAL(KIND=dp) :: TOL
4230
      INTEGER :: n
4231
      TYPE(Matrix_t) :: A
4232
!------------------------------------------------------------------------------
4233
      INTEGER, PARAMETER :: WORKN = 128
4234
      INTEGER :: i,j,k,l,istat, RowMin, RowMax
4235
      REAL(KIND=dp) :: NORMA
4236
      REAL(KIND=dp), POINTER CONTIG :: Values(:), ILUValues(:), CWork(:)
4237
      INTEGER, POINTER CONTIG :: Cols(:), Rows(:), Diag(:), &
4238
           ILUCols(:), ILURows(:), ILUDiag(:), IWork(:)
4239
      LOGICAL :: C(n)
4240
      REAL(KIND=dp) :: S(n), cptime, ttime, t
4241
!------------------------------------------------------------------------------
4242

4243
      ttime  = CPUTime()
4244
      cptime = 0.0d0
4245

4246
      Diag => A % Diag
4247
      Rows => A % Rows
4248
      Cols => A % Cols
4249
      IF(ASSOCIATED(A % PrecValues)) THEN
4250
        Values => A % PrecValues
4251
      ELSE
4252
        Values => A % Values
4253
      END IF
4254

4255
      ALLOCATE( A % ILURows(N+1),A % ILUDiag(N),STAT=istat )
4256
      IF ( istat /= 0 ) THEN
4257
         CALL Fatal( 'CRS_ILUT', 'Memory allocation error.' )
4258
      END IF
4259

4260
      ILURows => A % ILURows
4261
      ILUDiag => A % ILUDiag
4262

4263
      ALLOCATE( ILUCols( WORKN*N ),  ILUValues( WORKN*N ), STAT=istat )
4264
      IF ( istat /= 0 ) THEN
4265
         CALL Fatal( 'CRS_ILUT', 'Memory allocation error.' )
4266
      END IF
4267
!
4268
!     The factorization row by row:
4269
!     -----------------------------
4270
      ILURows(1) = 1
4271
      S =  0.0d0
4272
      C = .FALSE.
4273

4274
      DO i=1,n
4275
!
4276
!        Convert the current row to full form for speed,
4277
!        only flagging the nonzero entries:
4278
!        -----------------------------------------------
4279
         DO k=Rows(i), Rows(i+1) - 1
4280
            C(Cols(k)) = .TRUE.
4281
            S(Cols(k)) = Values(k)
4282
         END DO
4283
!
4284
!        Check bandwidth for speed, bandwidth optimization
4285
!        helps here ALOT, use it!
4286
!        -------------------------------------------------
4287
         RowMin = Cols(Rows(i))
4288
         RowMax = Cols(Rows(i+1)-1)
4289
!
4290
!        Here is the factorization part for the current row:
4291
!        ---------------------------------------------------
4292
         DO k=RowMin,i-1
4293
            IF ( C(k) ) THEN
4294
               IF ( ABS(ILUValues(ILUDiag(k))) > AEPS ) &
4295
                 S(k) = S(k) / ILUValues(ILUDiag(k)) 
4296
              
4297
               DO l=ILUDiag(k)+1, ILURows(k+1)-1
4298
                  j = ILUCols(l)
4299
                  IF ( .NOT. C(j) ) THEN
4300
                     C(j) = .TRUE.
4301
                     RowMax = MAX( RowMax,j )
4302
                  END IF
4303
                  S(j) = S(j) - S(k) * ILUValues(l)
4304
               END DO
4305
            END IF
4306
         END DO
4307
!
4308
!        This is the ILUT part, drop element ILU(i,j), if
4309
!        ABS(ILU(i,j)) <= NORM(A(i,:))*TOL:
4310
!        -------------------------------------------------
4311
         NORMA = SQRT( SUM( ABS(Values(Rows(i):Rows(i+1)-1))**2 ) )
4312

4313
         j = ILURows(i)-1
4314
         DO k=RowMin, RowMax
4315
            IF ( C(k) ) THEN
4316
               IF ( ABS(S(k)) >= TOL*NORMA .OR. k==i ) THEN
4317
                  j = j + 1
4318
                  ILUCols(j)   = k
4319
                  ILUValues(j) = S(k)
4320
                  IF ( k == i ) ILUDiag(i) = j
4321
               END IF
4322
               S(k) =  0.0d0
4323
               C(k) = .FALSE.
4324
            END IF
4325
         END DO
4326
         ILURows(i+1) = j + 1
4327
!
4328
!        Preparations for the next row:
4329
!        ------------------------------
4330
         IF ( i < N ) THEN
4331
!
4332
!           Check if still enough workspace:
4333
!           --------------------------------
4334
            IF ( SIZE(ILUCols) < ILURows(i+1) + N ) THEN
4335

4336
               t = CPUTime()
4337
!              k = ILURows(i+1) + MIN( WORKN, n-i ) * n
4338
               k = ILURows(i+1) + MIN( 0.75d0*ILURows(i+1), (n-i)*(1.0d0*n) )
4339
               ALLOCATE( IWork(k), STAT=istat )
4340
               IF ( istat /= 0 ) THEN
4341
                  CALL Fatal( 'CRS_ILUT', 'Memory allocation error.' )
4342
               END IF
4343
               DO j=1,ILURows(i+1)-1
4344
                 IWork(j) = ILUCols(j)
4345
               END DO
4346
               DEALLOCATE( ILUCols )
4347

4348
               ALLOCATE( CWork(k), STAT=istat )
4349
               IF ( istat /= 0 ) THEN
4350
                  CALL Fatal( 'CRS_ILUT', 'Memory allocation error.' )
4351
               END IF
4352
               DO j=1,ILURows(i+1)-1
4353
                 CWork(j) = ILUValues(j)
4354
               END DO
4355
               DEALLOCATE( ILUValues )
4356

4357
               ILUCols   => IWork
4358
               ILUValues => CWork
4359
               NULLIFY( IWork, CWork )
4360
               cptime = cptime + ( CPUTime() - t )
4361
!              PRINT*,'tot: ', CPUTime()-ttime, 'copy: ', cptime
4362
            END IF
4363
         END IF
4364
      END DO
4365
!
4366
!     Prescale the diagonal for the LU solve:
4367
!     ---------------------------------------
4368
      DO i=1,n
4369
         IF ( ABS(ILUValues(ILUDiag(i))) < AEPS ) THEN
4370
            ILUValues(ILUDiag(i)) = 1.0d0
4371
         ELSE
4372
            ILUValues(ILUDiag(i)) = 1.0d0 / ILUValues(ILUDiag(i))
4373
         END IF
4374
      END DO
4375

4376
      A % ILUCols   => ILUCols
4377
      A % ILUValues => ILUValues
4378
!------------------------------------------------------------------------------
4379
    END SUBROUTINE ComputeILUT
4380
!------------------------------------------------------------------------------
4381
  END FUNCTION CRS_ILUT
4382
!------------------------------------------------------------------------------
4383

4384

4385

4386
!------------------------------------------------------------------------------
4387
!>    Buids an incomplete (ILUT) factorization for an iterative solver
4388
!>    preconditioner. Complex matrix version.
4389
!------------------------------------------------------------------------------
4390
  FUNCTION CRS_ComplexILUT(A,TOL) RESULT(Status)
4391
!------------------------------------------------------------------------------
4392
    TYPE(Matrix_t) :: A  !< Structure holding input matrix, will also hold the factorization on exit.
4393
    REAL(KIND=dp), INTENT(IN) :: TOL  !< Drop toleranece: if ILUT(i,j) <= NORM(A(i,:))*TOL the value is dropped.
4394
    LOGICAL :: Status    !< Whether or not the factorization succeeded.
4395
!------------------------------------------------------------------------------
4396
    INTEGER :: n
4397
    REAL(KIND=dp) :: t
4398
!------------------------------------------------------------------------------
4399

4400
    CALL Info( 'CRS_ComplexILUT', 'ILU(T) (Complex), Starting factorization: ', Level=5 )
4401
    t = CPUTime()
4402

4403
    n = A % NumberOfRows / 2
4404

4405
    IF ( ASSOCIATED( A % CILUValues ) ) THEN
4406
       DEALLOCATE( A % ILURows, A % ILUCols, A % ILUDiag, A % CILUValues )
4407
    END IF
4408
!
4409
!   ... and then to the point:
4410
!   --------------------------
4411
    CALL ComplexComputeILUT( A, n, TOL )
4412
 
4413
!------------------------------------------------------------------------------
4414
    
4415
    WRITE( Message, * ) 'ILU(T) (Complex), NOF nonzeros: ',A % ILURows(n+1)
4416
    CALL Info( 'CRS_ComplexILUT', Message, Level=6 )
4417
    WRITE( Message, * ) 'ILU(T) (Complex), filling (%): ', &
4418
         FLOOR(A % ILURows(n+1)*(400.0d0/A % Rows(2*n+1)))
4419
    CALL Info( 'CRS_ComplexILUT', Message, Level=6 )
4420
    WRITE(Message,'(A,F8.2)') 'ILU(T) (Complex), Factorization ready at (s): ', CPUTime()-t
4421
    CALL Info( 'CRS_ComplexILUT', Message, Level=6 )
4422

4423
    Status = .TRUE.
4424
!------------------------------------------------------------------------------
4425

4426
  CONTAINS
4427

4428
!------------------------------------------------------------------------------
4429
    SUBROUTINE ComplexComputeILUT( A,n,TOL )
4430
!------------------------------------------------------------------------------
4431
      REAL(KIND=dp) :: TOL
4432
      INTEGER :: n
4433
      TYPE(Matrix_t) :: A
4434
!------------------------------------------------------------------------------
4435
      INTEGER, PARAMETER :: WORKN = 128
4436

4437
      INTEGER :: i,j,k,l,istat,RowMin,RowMax
4438
      REAL(KIND=dp) :: NORMA
4439

4440
      REAL(KIND=dp), POINTER CONTIG :: Values(:)
4441
      COMPLEX(KIND=dp), POINTER CONTIG :: ILUValues(:), CWork(:)
4442

4443
      INTEGER, POINTER CONTIG :: Cols(:), Rows(:), Diag(:), &
4444
           ILUCols(:), ILURows(:), ILUDiag(:), IWork(:)
4445

4446
      LOGICAL :: C(n)
4447
      COMPLEX(KIND=dp) :: S(n)
4448
!------------------------------------------------------------------------------
4449

4450
      Diag => A % Diag
4451
      Rows => A % Rows
4452
      Cols => A % Cols
4453
      IF (ASSOCIATED(A % PrecValues)) THEN
4454
        Values => A % PrecValues
4455
      ELSE
4456
        Values => A % Values
4457
      END IF
4458

4459
      ALLOCATE( A % ILURows(n+1),A % ILUDiag(n),STAT=istat )
4460
      IF ( istat /= 0 ) THEN
4461
         CALL Fatal( 'CRS_ComplexILUT', 'Memory allocation error.' )
4462
      END IF
4463

4464
      ILURows => A % ILURows
4465
      ILUDiag => A % ILUDiag
4466

4467
      ALLOCATE( ILUCols( WORKN*N ),  ILUValues( WORKN*N ), STAT=istat )
4468
      IF ( istat /= 0 ) THEN
4469
         CALL Fatal( 'CRS_ComplexILUT', 'Memory allocation error.' )
4470
      END IF
4471
!
4472
!     The factorization row by row:
4473
!     -----------------------------
4474
      ILURows(1) = 1
4475
      C = .FALSE.
4476
      S = CMPLX( 0.0d0, 0.0d0, KIND=dp )
4477

4478
      DO i=1,n
4479
!
4480
!        Convert the current row to full form for speed,
4481
!        only flagging the nonzero entries:
4482
!        -----------------------------------------------
4483
         DO k=Rows(2*i-1), Rows(2*i)-1,2
4484
            C((Cols(k)+1) / 2) = .TRUE.
4485
            S((Cols(k)+1) / 2) = CMPLX( Values(k), -Values(k+1), KIND=dp )
4486
         END DO
4487
!
4488
!        Check bandwidth for speed, bandwidth optimization
4489
!        helps here ALOT, use it!
4490
!        -------------------------------------------------
4491
         RowMin = (Cols(Rows(2*i-1)) + 1) / 2
4492
         RowMax = (Cols(Rows(2*i)-1) + 1) / 2
4493
!
4494
!        Here is the factorization part for the current row:
4495
!        ---------------------------------------------------
4496
         DO k=RowMin,i-1
4497
            IF ( C(k) ) THEN
4498
               IF ( ABS(ILUValues(ILUDiag(k))) > AEPS ) &
4499
                 S(k) = S(k) / ILUValues(ILUDiag(k)) 
4500
              
4501
               DO l=ILUDiag(k)+1, ILURows(k+1)-1
4502
                  j = ILUCols(l)
4503
                  IF ( .NOT. C(j) ) THEN
4504
                     C(j) = .TRUE.
4505
                     RowMax = MAX( RowMax,j )
4506
                  END IF
4507
                  S(j) = S(j) - S(k) * ILUValues(l)
4508
               END DO
4509
            END IF
4510
         END DO
4511

4512
!
4513
!        This is the ILUT part, drop element ILUT(i,j), if
4514
!        ABS(ILUT(i,j)) <= NORM(A(i,:))*TOL:
4515
!        -------------------------------------------------
4516
         NORMA = 0.0d0
4517
         DO k = Rows(2*i-1), Rows(2*i)-1, 2
4518
            NORMA = NORMA + Values(k)**2 + Values(k+1)**2
4519
         END DO
4520
         NORMA = SQRT(NORMA)
4521

4522
         j = ILURows(i)-1
4523
         DO k=RowMin, RowMax
4524
            IF ( C(k) ) THEN
4525
               IF ( ABS(S(k)) > TOL*NORMA .OR. k==i ) THEN
4526
                  j = j + 1
4527
                  ILUCols(j)   = k
4528
                  ILUValues(j) = S(k)
4529
                  IF ( k == i ) ILUDiag(i) = j
4530
               END IF
4531
               C(k) = .FALSE.
4532
               S(k) = CMPLX( 0.0d0, 0.0d0, KIND=dp )
4533
            END IF
4534
         END DO
4535
         ILURows(i+1) = j + 1
4536
!
4537
!        Preparations for the next row:
4538
!        ------------------------------
4539
         IF ( i < N ) THEN
4540
!
4541
!           Check if still enough workspace:
4542
!           --------------------------------
4543
            IF ( SIZE(ILUCols) < ILURows(i+1) + n ) THEN
4544
!              k = ILURows(i+1) + MIN( WORKN, n-i ) * n
4545
               k = ILURows(i+1) + MIN( 0.75d0*ILURows(i+1), (n-i)*(1.0d0*n) )
4546

4547
               ALLOCATE( IWork(k), STAT=istat )
4548
               IF ( istat /= 0 ) THEN
4549
                  CALL Fatal( 'CRS_ComplexILUT', 'Memory allocation error.' )
4550
               END IF
4551
               DO j=1,ILURows(i+1)-1
4552
                 IWork(j) = ILUCols(j)
4553
               END DO
4554
               DEALLOCATE( ILUCols )
4555

4556
               ALLOCATE( CWork(k), STAT=istat )
4557
               IF ( istat /= 0 ) THEN
4558
                  CALL Fatal( 'CRS_ComplexILUT', 'Memory allocation error.' )
4559
               END IF
4560
               DO j=1,ILURows(i+1)-1
4561
                 CWork(j) = ILUValues(j)
4562
               END DO
4563
               DEALLOCATE( ILUValues )
4564

4565
               ILUCols   => IWork
4566
               ILUValues => CWork
4567
            END IF
4568
         END IF
4569
      END DO
4570
!
4571
!     Prescale the diagonal for the LU solve:
4572
!     ---------------------------------------
4573
      DO i=1,n
4574
         IF ( ABS(ILUValues(ILUDiag(i))) < AEPS ) THEN
4575
            ILUValues(ILUDiag(i)) = 1.0d0
4576
         ELSE
4577
            ILUValues(ILUDiag(i)) = 1.0d0 / ILUValues(ILUDiag(i))
4578
         END IF
4579
      END DO
4580

4581
      A % ILUCols    => ILUCols
4582
      A % CILUValues => ILUValues
4583
      NULLIFY( ILUCols, ILUValues )
4584
!------------------------------------------------------------------------------
4585
    END SUBROUTINE ComplexComputeILUT
4586
!------------------------------------------------------------------------------
4587
  END FUNCTION CRS_ComplexILUT
4588
!------------------------------------------------------------------------------
4589

4590

4591

4592
!------------------------------------------------------------------------------
4593
!>    Incomplete factorization preconditioner solver for a CRS format matrix.
4594
!>    Matrix is accessed from a global variable GlobalMatrix. The real version.
4595
!------------------------------------------------------------------------------
4596
  SUBROUTINE CRS_LUPrecondition( u,v,ipar )
4597
!------------------------------------------------------------------------------
4598
    INTEGER, DIMENSION(*), INTENT(IN) :: ipar  !< structure holding info from (HUTIter-iterative solver package)
4599
    REAL(KIND=dp), DIMENSION(HUTI_NDIM), INTENT(IN) :: v   !< Right-hand-side vector
4600
    REAL(KIND=dp), DIMENSION(HUTI_NDIM), INTENT(OUT) :: u   !< Solution vector
4601
    
4602
    INTEGER :: i
4603

4604
    !$OMP PARALLEL DO
4605
    DO i=1,HUTI_NDIM
4606
       u(i) = v(i)
4607
    END DO
4608
    !$OMP END PARALLEL DO
4609
    CALL CRS_LUSolve( HUTI_NDIM,GlobalMatrix,u )
4610
  END SUBROUTINE CRS_LUPrecondition
4611
!------------------------------------------------------------------------------
4612

4613

4614

4615
!------------------------------------------------------------------------------
4616
!>    Incomplete factorization preconditioner solver for a CRS format matrix.
4617
!>    Matrix is accessed from a global variable GlobalMatrix. The Complex version.
4618
!------------------------------------------------------------------------------
4619
  SUBROUTINE CRS_ComplexLUPrecondition( u,v,ipar )
4620
!------------------------------------------------------------------------------
4621
    INTEGER, DIMENSION(*), INTENT(IN) :: ipar  !< structure holding info from (HUTIter-iterative solver package)
4622
    COMPLEX(KIND=dp), DIMENSION(HUTI_NDIM), INTENT(IN) :: v !< Right-hand-side vector
4623
    COMPLEX(KIND=dp), DIMENSION(HUTI_NDIM), INTENT(OUT) :: u !< Solution vector
4624

4625
    u = v
4626
    CALL CRS_ComplexLUSolve( HUTI_NDIM,GlobalMatrix,u )
4627
  END SUBROUTINE CRS_ComplexLUPrecondition
4628
!------------------------------------------------------------------------------
4629

4630

4631

4632
!------------------------------------------------------------------------------
4633
!>    Solve a system (Ax=b) after factorization A=LUD has been done. This
4634
!>    routine is meant as a part of  a preconditioner for an iterative solver.
4635
!------------------------------------------------------------------------------
4636
  SUBROUTINE CRS_LUSolve( N,A,b )
4637
!------------------------------------------------------------------------------
4638
    TYPE(Matrix_t), INTENT(IN) :: A  !< Structure holding input matrix
4639
    INTEGER, INTENT(IN) :: N   !< Size of the system
4640
    DOUBLE PRECISION :: b(n)   !< on entry the RHS vector, on exit the solution vector.
4641
!------------------------------------------------------------------------------
4642
    INTEGER :: i,j,k,l,row,col,nn
4643
    DOUBLE PRECISION :: s
4644
    DOUBLE PRECISION, POINTER CONTIG :: Values(:)
4645
    INTEGER, POINTER CONTIG :: Cols(:),Rows(:),Diag(:)
4646
!------------------------------------------------------------------------------
4647

4648
    Diag => A % ILUDiag
4649
    Rows => A % ILURows
4650
    Cols => A % ILUCols
4651
    Values => A % ILUValues
4652

4653
!
4654
!   if no ilu provided do diagonal solve:
4655
!   -------------------------------------
4656
    IF ( .NOT. ASSOCIATED( Values ) ) THEN
4657
       DO i=1,A % NumberOfRows
4658
         s = A % Values(A % Diag(i))
4659
         IF(s /= 0 ) b(i) = b(i) / s
4660
       END DO
4661
       RETURN
4662
    END IF
4663

4664
    IF ( A % Cholesky ) THEN
4665
      !
4666
      ! Forward substitute (solve z from Lz = b)
4667
      DO i=1,n
4668
        s = b(i)
4669
!DIR$ IVDEP
4670
        DO j=Rows(i),Diag(i)-1
4671
           s = s - Values(j) * b(Cols(j))
4672
        END DO
4673
        b(i) = s * Values(Diag(i))
4674
      END DO
4675

4676
      !
4677
      ! Backward substitute (solve x from L^Tx = z)
4678
      DO i=n,1,-1
4679
        b(i) = b(i) * Values(Diag(i))
4680
!DIR$ IVDEP
4681
        DO j=Rows(i),Diag(i)-1
4682
           b(Cols(j)) = b(Cols(j)) - Values(j) * b(i)
4683
        END DO
4684
      END DO
4685
    ELSE
4686
      !
4687
      ! Forward substitute (solve z from Lz = b)
4688
      DO i=1,n
4689
         s = b(i)
4690
!DIR$ IVDEP
4691
         DO j=Rows(i),Diag(i)-1
4692
            s = s - Values(j) * b(Cols(j))
4693
         END DO
4694
         b(i) = s 
4695
      END DO
4696

4697
      !
4698
      ! Backward substitute (solve x from UDx = z)
4699
      DO i=n,1,-1
4700
         s = b(i)
4701
!DIR$ IVDEP
4702
         DO j=Diag(i)+1,Rows(i+1)-1
4703
            s = s - Values(j) * b(Cols(j))
4704
         END DO
4705
         b(i) = Values(Diag(i)) * s
4706
      END DO
4707
    END IF
4708

4709
  END SUBROUTINE CRS_LUSolve
4710
!------------------------------------------------------------------------------
4711

4712

4713

4714
!------------------------------------------------------------------------------
4715
!>    Solve a complex system Ax=b after factorization A=LUD has been
4716
!>    done. This routine is meant as a part of  a preconditioner for an
4717
!>    iterative solver.
4718
!------------------------------------------------------------------------------
4719
  SUBROUTINE CRS_ComplexLUSolve( N,A,b )
4720
!------------------------------------------------------------------------------ 
4721
    TYPE(Matrix_t), INTENT(IN) :: A   !< Structure holding input matrix
4722
    INTEGER, INTENT(IN) :: N          !< Size of the system
4723
    COMPLEX(KIND=dp) :: b(N)          !< on entry the RHS vector, on exit the solution vector.
4724
!------------------------------------------------------------------------------
4725
    COMPLEX(KIND=dp), POINTER :: Values(:)
4726
    INTEGER :: i,j
4727
    COMPLEX(KIND=dp) :: x, s
4728
    INTEGER, POINTER :: Cols(:),Rows(:),Diag(:)
4729
!------------------------------------------------------------------------------
4730

4731
    Diag => A % ILUDiag
4732
    Rows => A % ILURows
4733
    Cols => A % ILUCols
4734
    Values => A % CILUValues
4735

4736
!
4737
!   if no ilu provided do diagonal solve:
4738
!   -------------------------------------
4739
    IF ( .NOT. ASSOCIATED( Values ) ) RETURN
4740

4741
    IF ( A % Cholesky ) THEN
4742
      !
4743
      ! Forward substitute
4744
      DO i=1,n
4745
         s = b(i)
4746
         DO j=Rows(i),Diag(i)-1
4747
            s = s - Values(j) * b(Cols(j))
4748
         END DO
4749
         b(i) = s * Values(Diag(i))
4750
      END DO
4751

4752
      !
4753
      ! Backward substitute
4754
      DO i=n,1,-1
4755
         b(i) = b(i) * Values(Diag(i))
4756
         DO j=Rows(i),Diag(i)-1
4757
           b(Cols(j)) = b(Cols(j)) - Values(j) * b(i)
4758
         END DO
4759
      END DO
4760
    ELSE
4761
      !
4762
      ! Forward substitute
4763
      DO i=1,n
4764
         s = b(i)
4765
         DO j=Rows(i),Diag(i)-1
4766
            s = s - Values(j) * b(Cols(j))
4767
         END DO
4768
         b(i) = s
4769
      END DO
4770

4771
      !
4772
      ! Backward substitute
4773
      DO i=n,1,-1
4774
         s = b(i)
4775
         DO j=Diag(i)+1,Rows(i+1)-1
4776
            s = s - Values(j) * b(Cols(j))
4777
         END DO
4778
         b(i) = Values(Diag(i)) * s
4779
      END DO
4780
    END IF
4781

4782
  END SUBROUTINE CRS_ComplexLUSolve
4783
!------------------------------------------------------------------------------
4784

4785

4786
!------------------------------------------------------------------------------
4787
!>    Matrix vector product (v = Au) for a matrix given in CRS format. The
4788
!>    matrix is accessed from a global variable GlobalMatrix.
4789
!------------------------------------------------------------------------------
4790
  SUBROUTINE CRS_MatrixVectorProd( u,v,ipar )
4791
!------------------------------------------------------------------------------
4792
    INTEGER, DIMENSION(*), INTENT(IN) :: ipar  !< Structure holding info HUTIter-iterative solver package
4793
    REAL(KIND=dp), INTENT(IN) :: u(*)  !< vector to multiply u
4794
    REAL(KIND=dp) :: v(*)              !< result vector
4795

4796
!------------------------------------------------------------------------------
4797
    INTEGER, POINTER  CONTIG :: Cols(:),Rows(:)
4798
    REAL(KIND=dp), POINTER  CONTIG :: Values(:)
4799

4800
#ifdef HAVE_MKL
4801
    INTERFACE
4802
      SUBROUTINE mkl_dcsrgemv(transa, m, a, ia, ja, x, y)
4803
        USE Types
4804
        CHARACTER :: transa
4805
        INTEGER :: m
4806
        REAL(KIND=dp) :: a(*)
4807
        INTEGER :: ia(*), ja(*)
4808
        REAL(KIND=dp) :: x(*), y(*)
4809
      END SUBROUTINE mkl_dcsrgemv
4810
    END INTERFACE
4811
#endif
4812

4813
    INTEGER :: i,j,l,n,ndeg
4814
    REAL(KIND=dp) :: s,r1,r2,r3,r4,r5
4815
    
4816
!------------------------------------------------------------------------------
4817

4818
    n = GlobalMatrix % NumberOfRows
4819
    Rows   => GlobalMatrix % Rows
4820
    Cols   => GlobalMatrix % Cols
4821
    Values => GlobalMatrix % Values
4822
    ndeg = GlobalMatrix % ndeg
4823
    
4824
    IF  ( GlobalMatrix % MatVecSubr /= 0 ) THEN
4825
      CALL MatVecSubrExt(GlobalMatrix % MatVecSubr, &
4826
          GlobalMatrix % SpMV, n,Rows,Cols,Values,u,v,0)
4827
      RETURN
4828
    END IF
4829

4830
    IF ( HUTI_EXTOP_MATTYPE == HUTI_MAT_NOTTRPSED ) THEN
4831
#ifdef HAVE_MKL
4832
      CALL mkl_dcsrgemv('N', n, Values, Rows, Cols, u, v)
4833
#else
4834

4835
    ! There may be a small structured block in the CRS matrix that is due to the problem
4836
    ! being initially vector valued. For example, in 3D elasticity we usually have dofs related
4837
    ! to (x,y,z) displacements following each other. Using this small dense block we may reduce
4838
    ! indirect memory addressing a little.
4839
    !-------------------------------------------------------------------------------------------
4840
      SELECT CASE( ndeg )
4841
      
4842
      CASE( 5, 10 )
4843
        !$omp parallel do private(j,l,r1,r2,r3,r4,r5)
4844
        DO i=1,n
4845
          r1 = 0.0_dp; r2 = 0.0_dp; r3 = 0.0_dp; r4 = 0.0_dp; r5 = 0.0_dp
4846
          !DIR$ IVDEP
4847
          DO j=Rows(i),Rows(i+1)-1,5
4848
            l = Cols(j)
4849
            r1 = r1 + u(l) * Values(j)
4850
            r2 = r2 + u(l+1) * Values(j+1)
4851
            r3 = r3 + u(l+2) * Values(j+2)
4852
            r4 = r4 + u(l+3) * Values(j+3)
4853
            r5 = r5 + u(l+4) * Values(j+4)
4854
          END DO
4855
          v(i) = r1 + r2 + r3 + r4 + r5
4856
        END DO
4857
        !$omp end parallel do
4858
        
4859
      CASE( 4, 8 )
4860
        !$omp parallel do private(j,l,r1,r2,r3,r4)
4861
        DO i=1,n
4862
          r1 = 0.0_dp; r2 = 0.0_dp; r3 = 0.0_dp; r4 = 0.0_dp
4863
          !DIR$ IVDEP
4864
          DO j=Rows(i),Rows(i+1)-1,4
4865
            l = Cols(j)
4866
            r1 = r1 + u(l) * Values(j)
4867
            r2 = r2 + u(l+1) * Values(j+1)
4868
            r3 = r3 + u(l+2) * Values(j+2)
4869
            r4 = r4 + u(l+3) * Values(j+3)
4870
          END DO
4871
          v(i) = r1 + r2 + r3 + r4
4872
        END DO
4873
        !$omp end parallel do
4874
        
4875
      CASE( 3, 6 )
4876
        !$omp parallel do shared(n,rows,cols,values) private(i,j,l,r1,r2,r3)
4877
        DO i=1,n
4878
          r1 = 0.0_dp; r2 = 0.0_dp; r3 = 0.0_dp
4879
          DO j=Rows(i),Rows(i+1)-1,3
4880
            l = Cols(j)
4881
            r1 = r1 + u(l) * Values(j)
4882
            r2 = r2 + u(l+1) * Values(j+1)
4883
            r3 = r3 + u(l+2) * Values(j+2)
4884
          END DO
4885
          v(i) = r1 + r2 + r3
4886
        END DO
4887
        !$omp end parallel do
4888
        
4889
      CASE( 2 )
4890
        !$omp parallel do private(j,l,r1,r2)
4891
        DO i=1,n
4892
          r1 = 0.0_dp; r2 = 0.0_dp
4893
          !DIR$ IVDEP
4894
          DO j=Rows(i),Rows(i+1)-1,2
4895
            l = Cols(j)
4896
            r1 = r1 + u(l) * Values(j)
4897
            r2 = r2 + u(l+1) * Values(j+1)
4898
          END DO
4899
          v(i) = r1 + r2
4900
        END DO
4901
        !$omp end parallel do
4902
        
4903
      CASE DEFAULT      
4904
        !$omp parallel do private(j,r1)
4905
        DO i=1,n
4906
          r1 = 0.0_dp
4907
          !DIR$ IVDEP
4908
          DO j=Rows(i),Rows(i+1)-1
4909
            r1 = r1 + u(Cols(j)) * Values(j)
4910
          END DO
4911
          v(i) = r1
4912
        END DO        
4913
        !$omp end parallel do      
4914

4915
      END SELECT
4916
#endif
4917
    ELSE
4918
      v(1:n) = 0.0d0
4919
      DO i=1,n
4920
        s = u(i)
4921
        DO j=Rows(i),Rows(i+1)-1
4922
          v(Cols(j)) = v(Cols(j)) + s * Values(j)
4923
        END DO
4924
      END DO
4925
    END IF
4926

4927
!    IF ( ASSOCIATED( GlobalMatrix % EMatrix ) ) THEN
4928
!       n = GlobalMatrix % EMatrix % NumberOFRows
4929
!       Rows   => GlobalMatrix % EMatrix % Rows
4930
!       Cols   => GlobalMatrix % EMatrix % Cols
4931
!       Values => GlobalMatrix % EMatrix % Values
4932
!
4933
!       allocate( w(n) )
4934
!       DO i=1,n
4935
!          s = 0.0d0
4936
!          DO j=Rows(i),Rows(i+1)-1
4937
!             s = s + Values(j) * u(Cols(j))
4938
!          END DO
4939
!          w(i) = s
4940
!       END DO
4941
!
4942
!       DO i=1,n
4943
!          s = w(i)
4944
!          DO j=Rows(i),Rows(i+1)-1
4945
!             v(Cols(j)) = v(Cols(j)) + s * Values(j)
4946
!          END DO
4947
!       END DO
4948
!       deallocate( w )
4949
!    END IF
4950

4951
  END SUBROUTINE CRS_MatrixVectorProd
4952
!------------------------------------------------------------------------------
4953

4954

4955
!------------------------------------------------------------------------------
4956
!>    Complex matrix vector product (v = Au) for a matrix given in
4957
!>    CRS format. The matrix is accessed from a global variable
4958
!>    GlobalMatrix.
4959
!------------------------------------------------------------------------------
4960
  SUBROUTINE CRS_ComplexMatrixVectorProd( u,v,ipar )
4961
!------------------------------------------------------------------------------
4962

4963
    INTEGER, DIMENSION(*), INTENT(IN) :: ipar      !< Structure holding info HUTIter-iterative solver package
4964
    COMPLEX(KIND=dp), INTENT(IN) :: u(HUTI_NDIM)      !< vector to multiply u
4965
    COMPLEX(KIND=dp) :: v(HUTI_NDIM)  !< result vector
4966

4967
!------------------------------------------------------------------------------
4968
    INTEGER, POINTER :: Cols(:),Rows(:)
4969
    INTEGER :: i,j,n
4970
    COMPLEX(KIND=dp) :: s,rsum
4971
    REAL(KIND=dp), POINTER :: Values(:)
4972
!------------------------------------------------------------------------------
4973

4974
    n = HUTI_NDIM
4975
    Rows   => GlobalMatrix % Rows
4976
    Cols   => GlobalMatrix % Cols
4977
    Values => GlobalMatrix % Values
4978

4979
    IF ( HUTI_EXTOP_MATTYPE == HUTI_MAT_NOTTRPSED ) THEN
4980
!$omp parallel do private(rsum,j,s)
4981
       DO i=1,n
4982
          rsum = CMPLX( 0.0d0, 0.0d0, KIND=dp )
4983
          DO j=Rows(2*i-1),Rows(2*i)-1,2
4984
             s = CMPLX( Values(j), -Values(j+1), KIND=dp )
4985
             rsum = rsum + s * u((Cols(j)+1)/2)
4986
          END DO
4987
          v(i) = rsum
4988
       END DO
4989
!$omp end parallel do
4990
    ELSE
4991
       v = CMPLX( 0.0d0, 0.0d0, KIND=dp )
4992
       DO i=1,n
4993
          rsum = u(i)
4994
!DIR$ IVDEP
4995
          DO j=Rows(2*i-1),Rows(2*i)-1,2
4996
             s = CMPLX( Values(j), -Values(j+1), KIND=dp )
4997
             v((Cols(j)+1)/2) = v((Cols(j)+1)/2) + s * rsum
4998
          END DO
4999
       END DO
5000
    END IF
5001

5002
  END SUBROUTINE CRS_ComplexMatrixVectorProd
5003
!------------------------------------------------------------------------------
5004

5005

5006
!------------------------------------------------------------------------------
5007
!> Check the matrix for correctness.
5008
!------------------------------------------------------------------------------
5009
  SUBROUTINE CRS_InspectMatrix( A )
5010
!------------------------------------------------------------------------------
5011
    TYPE(Matrix_t) :: A   !< Structure holding the matrix
5012
!------------------------------------------------------------------------------
5013

5014
    REAL(KIND=dp), POINTER :: Values(:)
5015
    INTEGER, POINTER :: Cols(:),Rows(:),Diag(:)
5016
    INTEGER :: i,j,k,k2,n,ColMin,ColMax,RowMin,RowMax,ColN,RowN,ValN,UnSym
5017
    REAL(KIND=dp) :: ValMin,ValMax,TotalSum,DiagSum
5018
    LOGICAL :: Hit
5019

5020
    Diag   => A % Diag
5021
    Rows   => A % Rows
5022
    Cols   => A % Cols
5023
    Values => A % Values
5024

5025
    n = A % NumberOfRows
5026
    RowMin = MINVAL( Rows )
5027
    RowMax = MAXVAL( Rows )
5028
    RowN = SIZE( Rows )
5029

5030
    PRINT *,'Rows (size '//I2S(RowN)//') range:',RowMin, RowMax
5031

5032
    IF( RowMin < 1 ) THEN
5033
      PRINT *,'Outliers:'
5034
      j =  0
5035
      DO i=1,RowN
5036
        IF( Rows(i) < 1 ) THEN
5037
          j = j + 1
5038
          IF( j < 10 ) PRINT *,'i',i,Rows(i)
5039
        END IF
5040
      END DO
5041
      IF( j > 0 ) THEN
5042
        PRINT *,'Number of row outliers:',j
5043
      END IF
5044
    END IF
5045

5046
    ColMin = MINVAL( Cols )
5047
    ColMax = MAXVAL( Cols )
5048
    ColN = SIZE( Cols )
5049
    PRINT *,'Cols (size '//I2S(ColN)//') range:',ColMin, ColMax
5050
    
5051
    IF( ColMin < 1 ) THEN
5052
      PRINT *,'Outliers:'
5053
      j = 0
5054
      DO i=1,ColN
5055
        IF( Cols(i) < 1 ) THEN
5056
          j = j + 1
5057
          IF( j < 10 )PRINT *,'i',i,Cols(i)
5058
        END IF
5059
      END DO
5060
      IF( j > 0 ) THEN
5061
        PRINT *,'Number of column outliers:',j
5062
      END IF
5063
    END IF
5064

5065
    ValMin = MINVAL( Values )
5066
    ValMax = MAXVAL( Values )
5067
    ValN = SIZE( Values )
5068
    TotalSum = SUM( Values )
5069

5070
    PRINT *,'Values (size '//I2S(ValN)//') range:',ValMin,ValMax,TotalSum
5071
    
5072
    IF( ColN /= RowMax - 1  ) THEN
5073
      PRINT *,'Conflicting max row index :',n,RowMax-1
5074
    END IF
5075
    IF( N /= ColMax  ) THEN
5076
      PRINT *,'Conflicting max row index:',N,ColMax
5077
    END IF
5078
    IF( ColN /= ValN  ) THEN
5079
      PRINT *,'Conflicting size of Values vs. Cols:',ValN,ColN
5080
    END IF
5081

5082
    PRINT *,'Checking Diag vector:'
5083
    j = 0
5084
    DO i=1,n
5085
      IF( i /= Cols(Diag(i)) ) THEN
5086
        j = j + 1
5087
        IF( j < 10 ) PRINT *,'diag:',i,Cols(Diag(i)) 
5088
      END IF
5089
    END DO
5090
    IF( j > 0 ) THEN
5091
      PRINT *,'Number of erroneous diag entries:',j
5092
    ELSE
5093
      DiagSum = SUM( Values(Diag) )
5094
      PRINT *,'diagonal sum:',DiagSum
5095

5096
! ratios are used to study cluster MG method, for example
5097
      PRINT *,'ratios:',TotalSum/DiagSum,(TotalSum-DiagSum)/DiagSum
5098
      
5099
      TotalSum = SUM( ABS( Values ) )
5100
      DiagSum = SUM( ABS( Values(Diag) ) )
5101
      PRINT *,'abs ratios:',TotalSum/DiagSum,(TotalSum-DiagSum)/DiagSum
5102

5103
    END IF
5104
    
5105
    UnSym = 0
5106
    DO i=1,n
5107
      DO k=Rows(i),Rows(i+1)-1
5108
        j = Cols(k)
5109
        Hit = .FALSE.
5110
        DO k2=Rows(j),Rows(j+1)-1
5111
          IF( Cols(k2) == i ) THEN
5112
            Hit = .TRUE.
5113
            EXIT
5114
          END IF
5115
        END DO
5116
        IF( .NOT. Hit ) THEN
5117
          UnSym = UnSym + 1
5118
          IF( UnSym < 10 ) PRINT *,'unsym:',i,j
5119
        END IF
5120
      END DO
5121
    END DO
5122

5123
    PRINT *,'Number of unsymmetric entries in matrix topology:',UnSym
5124

5125

5126
  END SUBROUTINE CRS_InspectMatrix
5127

5128

5129
!------------------------------------------------------------------------------
5130
!> Eliminate redundant entries in matrix A by summing same (i,j) combinations
5131
!> together. It is assumed that the same entries are on the same CRS row.
5132
!------------------------------------------------------------------------------
5133
  SUBROUTINE CRS_PackMatrix( A )
5134
!------------------------------------------------------------------------------
5135
    TYPE(Matrix_t) :: A     !< Structure holding the matrix
5136
!------------------------------------------------------------------------------ 
5137
    INTEGER :: i,j,k,k2,n,rowi,nofs0
5138
    INTEGER, ALLOCATABLE :: LocalCols(:), ColIndex(:)
5139
    INTEGER, POINTER :: Cols(:),Rows(:),Diag(:)
5140
    REAL(KIND=dp), POINTER :: Values(:), TValues(:)
5141
    REAL(KIND=dp), ALLOCATABLE :: LocalValues(:), LocalTValues(:)
5142
!------------------------------------------------------------------------------
5143

5144
    IF(A % NumberOfRows == 0 ) RETURN
5145

5146
    Rows    => A % Rows
5147
    Cols    => A % Cols
5148
    Diag    => A % Diag
5149
    Values  => A % Values
5150
    TValues => A % TValues
5151

5152
    nofs0 = Rows(A % NumberOfRows+1)-1
5153

5154
    n = 0
5155
    DO i=1,A % NumberOfRows
5156
      n = MAX( n, Rows(i+1)+1-Rows(i) )
5157
    END DO
5158

5159
    ALLOCATE( LocalCols(n), LocalValues(n), LocalTValues(n) )
5160
    LocalCols = 0
5161
    LocalValues = 0.0_dp
5162

5163
    ALLOCATE(ColIndex(MAXVAL(Cols)))
5164
    ColIndex = 0
5165

5166
    k2 = 0
5167
    DO i=1,A % NumberOfRows
5168
      Rowi = k2+1
5169

5170
      ! Memorize the matrix row
5171
      DO k=1,Rows(i+1)-Rows(i)
5172
        LocalCols(k)   = Cols(Rows(i)+k-1)
5173
        LocalValues(k) = Values(Rows(i)+k-1)
5174
        IF(ASSOCIATED(TValues)) &
5175
          LocalTValues(k) = TValues(Rows(i)+k-1)
5176
      END DO
5177

5178
      ! Pack the matrix row 
5179
      DO k=1,Rows(i+1)-Rows(i)
5180
        j = LocalCols(k) 
5181
        IF( ColIndex(j) == 0 ) THEN
5182
          k2 = k2 + 1
5183
          ColIndex(j) = k2
5184
          Cols(k2) = Cols(Rows(i)+k-1)
5185
          Values(k2) = LocalValues(k)
5186
          IF(ASSOCIATED(TValues)) &
5187
            TValues(k2) = LocalTValues(k)
5188
        ELSE
5189
          Values(ColIndex(j)) = Values(ColIndex(j)) + LocalValues(k)
5190
          IF(ASSOCIATED(TValues)) &
5191
            TValues(ColIndex(j)) = TValues(ColIndex(j)) + LocalTValues(k)
5192
        END IF
5193
      END DO
5194
      
5195
      ! Nullify the index table
5196
      DO k=1,Rows(i+1)-Rows(i)
5197
        j = LocalCols(k) 
5198
        ColIndex(j) = 0 
5199
      END DO
5200
      Rows(i) = Rowi
5201
    END DO
5202
    Rows(i) = k2+1
5203

5204
    CALL Info('CRS_PackMatrix','Number of summed-up matrix entries: '&
5205
        //I2S(nofs0-k2),Level=8)
5206

5207
  END SUBROUTINE CRS_PackMatrix
5208
!------------------------------------------------------------------------------
5209

5210

5211
!------------------------------------------------------------------------------
5212
!> At first call register the matrix topology.
5213
!> At second round change the matrix topology of the vectors.
5214
!> Without change in topology all matrix operations with BulkValues,
5215
!> MassValues, DampValues etc. would break down. 
5216
!------------------------------------------------------------------------------
5217
  SUBROUTINE CRS_ChangeTopology( A, Init )
5218
!------------------------------------------------------------------------------
5219
    TYPE(Matrix_t) :: A     !< Structure holding the matrix
5220
    LOGICAL :: Init
5221
!------------------------------------------------------------------------------ 
5222
    LOGICAL, SAVE :: InitDone = .FALSE.
5223
    INTEGER, ALLOCATABLE, SAVE :: Rows0(:), Cols0(:)
5224
    REAL(KIND=dp), POINTER :: Aold(:), Anew(:)
5225
    INTEGER, SAVE :: n0
5226
    INTEGER :: i,j,k,j2,k2,ivec,n
5227
    
5228

5229
    IF( A % NumberOfRows == 0 ) RETURN
5230
        
5231
    IF(Init) THEN
5232
      IF(InitDone) THEN
5233
        CALL Warn('CRS_ChangeTopology','We already have initialized Cols0 and Rows0!')
5234
        DEALLOCATE(Cols0,Rows0)
5235
      END IF
5236
      n0 = SIZE(A % Cols)
5237
      CALL Info('CRS_ChangeTopology','Original matrix non-zeros: '//I2S(n0),Level=12)
5238
      
5239
      ALLOCATE( Cols0(n0), Rows0( SIZE( A % Rows ) ) )
5240
      Cols0 = A % Cols
5241
      Rows0 = A % Rows
5242
      InitDone = .TRUE.
5243
    ELSE
5244
      n = SIZE( A % Cols )
5245
      
5246
      IF( n == n0 ) THEN
5247
        IF( ALL( Cols0 == A % Cols ) ) THEN
5248
          CALL Info('CRS_ChangeTopology','Topology is unaltered!',Level=20)
5249
          DEALLOCATE(Cols0,Rows0)
5250
          InitDone = .FALSE.
5251
          RETURN
5252
        END IF
5253
      END IF
5254

5255
      IF( SIZE(A % Rows) /= SIZE(Rows0) ) THEN
5256
        CALL Fatal('CRS_ChangeTopology','This routine assumes constant number of rows!')
5257
      END IF
5258
      
5259
      CALL Info('CRS_ChangeTopology','New matrix non-zeros: '//I2S(n),Level=12)
5260
      
5261
      DO ivec=1,5
5262
        NULLIFY(Aold)
5263
        SELECT CASE(ivec)
5264
        CASE( 1 )
5265
          Aold => A % BulkValues
5266
        CASE( 2 )
5267
          Aold => A % MassValues
5268
        CASE( 3 )
5269
          Aold => A % DampValues
5270
        CASE( 4 )
5271
          Aold => A % BulkMassValues
5272
        CASE( 5 )
5273
          Aold => A % BulkDampValues
5274
        END SELECT
5275
        
5276
        IF( .NOT. ASSOCIATED(Aold) ) CYCLE
5277
        
5278
        NULLIFY(Anew)
5279
        ALLOCATE(Anew(n))
5280
        Anew = 0.0_dp
5281
        
5282
        DO i=1,A % NumberOfRows
5283
          DO j = Rows0(i), Rows0(i+1)-1
5284
            k = Cols0(j) 
5285
            DO j2 = A % Rows(i), A % Rows(i+1)-1
5286
              k2 = A % Cols(j2)
5287
              IF( k == k2 ) THEN
5288
                Anew(j2) = Aold(j)
5289
                EXIT
5290
              END IF
5291
            END DO
5292
          END DO
5293
        END DO
5294
        
5295
        DEALLOCATE( Aold ) 
5296

5297
        SELECT CASE(ivec)
5298
        CASE( 1 )
5299
          A % BulkValues => Anew
5300
        CASE( 2 )
5301
          A % MassValues => Anew
5302
        CASE( 3 )
5303
          A % DampValues => Anew 
5304
        CASE( 4 )
5305
          A % BulkMassValues => Anew
5306
        CASE( 5 )
5307
          A % BulkDampValues => Anew 
5308
        END SELECT
5309

5310
        CALL Info('CRS_ChangeTopology','Done changing matrix '//I2S(n)//' topology',Level=20)
5311

5312
        
5313
      END DO
5314
            
5315
      DEALLOCATE(Cols0,Rows0)
5316
      InitDone = .FALSE.
5317
      
5318
      A % ndeg = -1 
5319
      CALL Info('CRS_ChangeTopology','Matrix topology changed',Level=30)
5320
    END IF
5321
          
5322
  END SUBROUTINE CRS_ChangeTopology
5323

5324
!------------------------------------------------------------------------------
5325
!> Check that matrix has a repeating block of size "dofs" that can be
5326
!> utilized in Matrix-Vector products, for example. 
5327
!------------------------------------------------------------------------------
5328
  FUNCTION CRS_CheckStructuredDofs( A, dofs) RESULT ( Failed )
5329
!------------------------------------------------------------------------------
5330
    TYPE(Matrix_t), INTENT(IN) :: A                !< Structure holding matrix
5331
    INTEGER :: dofs                                !< Size of dense block to be tested
5332
    LOGICAL :: Failed
5333
!------------------------------------------------------------------------------
5334
    INTEGER, POINTER  CONTIG :: Cols(:),Rows(:)
5335
    INTEGER :: i,j,k,n,m
5336
!------------------------------------------------------------------------------
5337
    n = A % NumberOfRows
5338
    Rows   => A % Rows
5339
    Cols   => A % Cols
5340

5341
    Failed = .FALSE.
5342
    DO i=1,n
5343
      DO j=Rows(i),Rows(i+1)-1,dofs
5344
        DO k=1,dofs-1
5345
          IF( Cols(j+k)-Cols(j) /= k ) THEN
5346
            Failed = .TRUE.
5347
            EXIT
5348
          END IF
5349
        END DO
5350
      END DO
5351
    END DO
5352
      
5353
  END FUNCTION CRS_CheckStructuredDofs
5354

5355

5356

5357
  
5358
END MODULE CRSMatrix
5359
!------------------------------------------------------------------------------
5360

5361
!> \} ElmerLib
5362

5363
Product

Resources

Company
