1
2!
3!  Dalton, a molecular electronic structure program
4!  Copyright (C) by the authors of Dalton.
5!
6!  This program is free software; you can redistribute it and/or
7!  modify it under the terms of the GNU Lesser General Public
8!  License version 2.1 as published by the Free Software Foundation.
9!
10!  This program is distributed in the hope that it will be useful,
11!  but WITHOUT ANY WARRANTY; without even the implied warranty of
12!  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
13!  Lesser General Public License for more details.
14!
15!  If a copy of the GNU LGPL v2.1 was not distributed with this
16!  code, you can obtain one at https://www.gnu.org/licenses/old-licenses/lgpl-2.1.en.html.
17!
18!
19C
20*---------------------------------------------------------------------*
21c/* Deck CC_KMATRIX */
22*=====================================================================*
23      SUBROUTINE CC_KMATRIX(IKTRAN, NKTRAN, LISTA, LISTB, IOPTRES,
24     &                    FILKMA, IKDOTS, KCONS, MXVEC, WORK, LWORK)
25*---------------------------------------------------------------------*
26*
27*    Purpose: batched loop over K matrix transformations
28*             (needed if the number of transformations exceeds the
29*              limit MAXSIM defined on ccsdio.h )
30*
31*    CCMM JK+OC, modyfied CC_BMATRIX
32*=====================================================================*
33#if defined (IMPLICIT_NONE)
34      IMPLICIT NONE
35#else
36#  include "implicit.h"
37#endif
38#include "priunit.h"
39#include "maxorb.h"
40#include "ccsdio.h"
41
42      LOGICAL LOCDBG
43      PARAMETER (LOCDBG = .FALSE.)
44
45      CHARACTER*(*) LISTA, LISTB, FILKMA
46      INTEGER IOPTRES
47      INTEGER NKTRAN, MXVEC, LWORK
48      INTEGER IKTRAN(3,NKTRAN)
49      INTEGER IKDOTS(MXVEC,NKTRAN)
50
51#if defined (SYS_CRAY)
52      REAL WORK(LWORK)
53      REAL KCONS(MXVEC,NKTRAN)
54#else
55      DOUBLE PRECISION WORK(LWORK)
56      DOUBLE PRECISION KCONS(MXVEC,NKTRAN)
57#endif
58
59      INTEGER MAXKTRAN, NTRAN, ISTART, IBATCH, NBATCH
60
61      MAXKTRAN = MAXSIM
62
63      NBATCH = (NKTRAN+MAXKTRAN-1)/MAXKTRAN
64
65      IF (LOCDBG) THEN
66        WRITE (LUPRI,*) 'Batching over K matrix transformations:'
67        WRITE (LUPRI,*) 'nb. of batches needed:', NBATCH
68      END IF
69
70      DO IBATCH = 1, NBATCH
71        ISTART = (IBATCH-1) * MAXKTRAN + 1
72        NTRAN  = MIN(NKTRAN-(ISTART-1),MAXKTRAN)
73
74        IF (LOCDBG) THEN
75          WRITE (LUPRI,*) 'Batch No.:',IBATCH
76          WRITE (LUPRI,*) 'start at :',ISTART
77          WRITE (LUPRI,*) '# transf.:',NTRAN
78        END IF
79
80        CALL CC_KMAT(IKTRAN(1,ISTART), NTRAN,
81     &                LISTA, LISTB, IOPTRES, FILKMA,
82     &                IKDOTS(1,ISTART), KCONS(1,ISTART),
83     &                MXVEC, WORK,LWORK)
84
85      END DO
86
87      RETURN
88      END
89
90*---------------------------------------------------------------------*
91*              END OF SUBROUTINE CC_KMATRIX                           *
92*---------------------------------------------------------------------*
93
94*---------------------------------------------------------------------*
95c/* Deck CC_BMAT */
96*=====================================================================*
97      SUBROUTINE CC_KMAT(IBTRAN, NBTRAN, LISTA, LISTB, IOPTRES,
98     &                    FILBMA, IBDOTS, BCONS, MXVEC, WORK, LWORK )
99*---------------------------------------------------------------------*
100*
101*             The linear transformations are calculated for a list
102*             of bar{T^A} vectors and a list of bar{T^B} vectors:
103*
104*                LISTA       -- type of bar{T^A} vectors
105*                LISTB       -- type of bar{T^B} vectors
106*                IBTRAN(1,*) -- indeces of bar{T^A} vectors
107*                IBTRAN(2,*) -- indeces of bar{T^B} vectors
108*                IBTRAN(3,*) -- indeces or addresses of result vectors
109*                NBTRAN      -- number of requested transformations
110*                FILBMA      -- file name / list type of result vectors
111*                               or list type of vectors to be dotted on
112*                IBDOTS      -- indeces of vectors to be dotted on
113*                BCONS       -- contains the dot products on return
114*
115*    return of the result vectors:
116*
117*           IOPTRES = 0 :  all result vectors are written to a direct
118*                          access file, FILBMA is used as file name
119*                          the start addresses of the vectors are
120*                          returned in IBTRAN(3,*)
121*
122*           IOPTRES = 1 :  the vectors are kept and returned in WORK
123*                          if possible, start addresses returned in
124*                          IBTRAN(3,*). N.B.: if WORK is not large
125*                          enough IOPTRES is automatically reset to 0!!
126*
127*           IOPTRES = 3 :  each result vector is written to its own
128*                          file by a call to CC_WRRSP, FILBMA is used
129*                          as list type and IBTRAN(3,*) as list index
130*                          NOTE that IBTRAN(3,*) is in this case input!
131*
132*           IOPTRES = 4 :  each result vector is added to a vector on
133*                          file by a call to CC_WARSP, FIBBMA is used
134*                          as list type and IBTRAN(3,*) as list index
135*                          NOTE that IBTRAN(3,*) is in this case input!
136*
137*           IOPTRES = 5 :  the result vectors are dotted on a array
138*                          of vectors, the type of the arrays given
139*                          by FILBMA and the indeces from IBDOTS
140*                          the result of the dot products is returned
141*                          in the BCONS array
142*
143*
144*           CCMM JK+OC, modyfied version of CC_BMAT
145*
146*=====================================================================*
147      USE PELIB_INTERFACE, ONLY: USE_PELIB, PELIB_IFC_QRTRANSFORMER
148#if defined (IMPLICIT_NONE)
149      IMPLICIT NONE
150#else
151#  include "implicit.h"
152#endif
153#include "priunit.h"
154#include "ccsdinp.h"
155#include "ccsdsym.h"
156#include "ccsections.h"
157#include "maxorb.h"
158#include "mxcent.h"
159#include "ccsdio.h"
160#include "ccorb.h"
161#include "iratdef.h"
162#include "eribuf.h"
163#include "ccslvinf.h"
164#include "qm3.h"
165#include "second.h"
166
167* local parameters:
168      CHARACTER MSGDBG*(16)
169      PARAMETER (MSGDBG='[debug] CC_KMAT> ')
170
171      LOGICAL LOCDBG, LSAME
172      PARAMETER (LOCDBG = .FALSE.)
173
174      INTEGER KDUM
175      PARAMETER( KDUM = +99 999 999 ) ! dummy address for work space
176
177
178      INTEGER LUBMAT
179
180      CHARACTER*(*) LISTA, LISTB, FILBMA
181
182      INTEGER IOPTRES
183      INTEGER NBTRAN, MXVEC, LWORK
184      INTEGER IBTRAN(3,NBTRAN)
185      INTEGER IBDOTS(MXVEC,NBTRAN)
186
187#if defined (SYS_CRAY)
188      REAL WORK(LWORK)
189      REAL ZERO, ONE, TWO
190      REAL DUM, XNORM, DUMMY
191      REAL BCONS(MXVEC,NBTRAN)
192      REAL FACTSLV
193#else
194      DOUBLE PRECISION WORK(LWORK)
195      DOUBLE PRECISION ZERO, ONE, TWO
196      DOUBLE PRECISION DUM, XNORM, DUMMY
197      DOUBLE PRECISION BCONS(MXVEC,NBTRAN)
198      DOUBLE PRECISION FACTSLV
199      DOUBLE PRECISION TAL1, TAL2, RHO1N, RHO2N
200#endif
201      PARAMETER (ZERO = 0.0d0, ONE = 1.0d0, TWO = 2.0d0)
202
203      CHARACTER*(10) MODEL, MODELW, CDUMMY
204      CHARACTER*8 LABEL
205      CHARACTER RSPTYP*(1)
206      INTEGER NBATCH
207
208
209      INTEGER ITRAN, IDLSTA, IDLSTB, IOPT
210      INTEGER ISYMA, ISYMB, ISYMAB
211      INTEGER IBATCH,IADRTH
212      INTEGER KEND1, LEN, LENALL
213      INTEGER KEND2, LWRK2, KEND3, LWRK3, LWRK1
214      INTEGER IVEC
215      INTEGER KT2AMPA, KTHETA0
216      INTEGER KTHETA1, KTHETA2, KT1AMPA, KT1AMPB
217      INTEGER IOPTW, IDUMMY
218      INTEGER KTGB, KETA, KETA1, KETA2, NAMPF
219      INTEGER KTGA
220
221* external functions:
222      INTEGER ILSTSYM
223
224#if defined (SYS_CRAY)
225      REAL DDOT, DTIME, TIMALL, TIMTRN
226#else
227      DOUBLE PRECISION DDOT, DTIME, TIMALL, TIMTRN
228#endif
229
230*---------------------------------------------------------------------*
231* begin:
232*---------------------------------------------------------------------*
233      IF (LOCDBG) THEN
234        Call AROUND('ENTERED CC_KMAT')
235        WRITE (LUPRI,*) 'LISTA : ',LISTA
236        WRITE (LUPRI,*) 'LISTB : ',LISTB
237        WRITE (LUPRI,*) 'FILBMA: ',FILBMA
238        WRITE (LUPRI,*) 'NKTRAN: ',NBTRAN
239        WRITE (LUPRI,*) 'IOPTRES:',IOPTRES
240        CALL FLSHFO(LUPRI)
241      END IF
242
243      IF (CCSDT) THEN
244        WRITE(LUPRI,'(/1x,a)') 'K matrix transformations not '
245     &          //'implemented for triples yet...'
246        CALL QUIT('Triples not implemented for K matrix '//
247     &            'transformations')
248      END IF
249
250      IF ( .not. (CCS .or. CC2 .or. CCSD) ) THEN
251        WRITE(LUPRI,'(/1x,a)') 'CC_KMAT called for a Coupled Cluster '
252     &          //'method not implemented in CC_KMAT...'
253        CALL QUIT('Unknown CC method in CC_KMAT.')
254      END IF
255
256      IF (.NOT. DUMPCD) THEN
257        WRITE(LUPRI,*) 'DUMPCD = ',DUMPCD
258        WRITE(LUPRI,*) 'CC_KMAT requires DUMPCD=.TRUE.'
259        CALL QUIT('DUMPCD=.FALSE. , CC_KMAT requires DUMPCD=.TRUE.')
260      END IF
261
262      IF (.NOT. RSPIM) THEN
263        WRITE(LUPRI,*) 'RSPIM = ',RSPIM
264        WRITE(LUPRI,*) 'CC_BMAT requires RSPIM=.TRUE.'
265        CALL QUIT('RSPIM=.FALSE. , CC_KMAT requires RSPIM=.TRUE.')
266      END IF
267
268      IF (ISYMOP .NE. 1) THEN
269        WRITE(LUPRI,*) 'ISYMOP = ',ISYMOP
270        WRITE(LUPRI,*) 'CC_KMAT is not implemented for ISYMOP.NE.1'
271        CALL QUIT('CC_KMAT is not implemented for ISYMOP.NE.1')
272      END IF
273
274      IF (NBTRAN .GT. MAXSIM) THEN
275        WRITE(LUPRI,*) 'NBTRAN = ', NBTRAN
276        WRITE(LUPRI,*) 'MAXSIM = ', MAXSIM
277        WRITE(LUPRI,*) 'number of requested transformation is larger'
278        WRITE(LUPRI,*) 'than the maximum number of allowed ',
279     &                 'simultaneous transformation.'
280        WRITE(LUPRI,*) 'Error in CC_KMAT: NBTRAN is larger than MAXSIM.'
281        CALL QUIT('Error in CC_KMAT: NBTRAN is larger than MAXSIM.')
282      END IF
283
284      IF (IPRINT.GT.0) THEN
285
286         WRITE (LUPRI,'(//1X,A1,50("="),A1)')'+','+'
287
288         WRITE (LUPRI,'(1x,A52)')
289     &         '|        K MATRIX TRANSFORMATION SECTION           |'
290
291         IF (IOPTRES.EQ.3) THEN
292            WRITE (LUPRI,'(1X,A52)')
293     &         '|          (result is written to file)             |'
294         ELSE IF (IOPTRES.EQ.4) THEN
295            WRITE (LUPRI,'(1X,A52)')
296     &         '|     (result is added to a vector on file)        |'
297         ELSE IF (IOPTRES.EQ.5) THEN
298            WRITE (LUPRI,'(1X,A52)')
299     &         '|    (result used to calculate dot products)       |'
300         END IF
301
302         WRITE (LUPRI,'(1X,A1,50("-"),A1)') '+','+'
303
304      END IF
305
306* initialize timings:
307      TIMALL  = SECOND()
308
309* set option and model to write vectors to file:
310      IF (CCS) THEN
311         MODELW = 'CCS       '
312         IOPTW  = 1
313      ELSE IF (CC2) THEN
314         MODELW = 'CC2       '
315         IOPTW  = 3
316      ELSE IF (CCSD) THEN
317         MODELW = 'CCSD      '
318         IOPTW  = 3
319      ELSE
320         CALL QUIT('Unknown coupled cluster model in CC_KMAT.')
321      END IF
322
323
324* check return option for the result vectors:
325      LUBMAT = -1
326      IF (IOPTRES .EQ. 0 .OR. IOPTRES .EQ. 1) THEN
327         CALL WOPEN2(LUBMAT, FILBMA, 64, 0)
328      ELSE IF (IOPTRES .EQ. 3 .OR. IOPTRES .EQ. 4) THEN
329         CONTINUE
330      ELSE IF (IOPTRES .EQ. 5) THEN
331         IF (MXVEC*NBTRAN.NE.0) CALL DZERO(BCONS,MXVEC*NBTRAN)
332      ELSE
333         CALL QUIT('Illegal value of IOPTRES in CC_KMAT.')
334      END IF
335C
336C----------------------------------------------
337C     If all models are SPC
338C     -> RETURN from CCMM_TGB:
339C----------------------------------------------
340C
341      IF (LOSPC) RETURN
342C
343*=====================================================================*
344* calculate K matrix transformations:
345*=====================================================================*
346      IADRTH = 1
347      DO ITRAN = 1, NBTRAN
348
349        IDLSTA = IBTRAN(1,ITRAN)
350        IDLSTB = IBTRAN(2,ITRAN)
351
352        ISYMA  = ILSTSYM(LISTA,IDLSTA)
353        ISYMB  = ILSTSYM(LISTB,IDLSTB)
354        ISYMAB = MULD2H(ISYMA,ISYMB)
355
356        TIMTRN = SECOND()
357
358*---------------------------------------------------------------------*
359* allocate work space for the result vector:
360*---------------------------------------------------------------------*
361        IF (CCS) THEN
362          KTHETA1 = 1
363          KTHETA2 = KDUM
364          KEND1   = KTHETA1 + NT1AM(ISYMAB)
365          LWRK1 = LWORK - KEND1
366          CALL DZERO(WORK(KTHETA1),NT1AM(ISYMAB))
367        ELSE
368          KTHETA1 = 1
369          KTHETA2 = KTHETA1 + NT1AM(ISYMAB)
370          KEND1   = KTHETA2 + NT2AM(ISYMAB)
371          LWRK1 = LWORK - KEND1
372          CALL DZERO(WORK(KTHETA1),NT1AM(ISYMAB))
373          CALL DZERO(WORK(KTHETA2),NT2AM(ISYMAB))
374        END IF
375
376        IF (LOCDBG) THEN
377         WRITE (LUPRI,*) 'K matrix transformation for ITRAN,',ITRAN
378         WRITE (LUPRI,*) 'IADRTH:',IADRTH
379         WRITE (LUPRI,*) 'LISTA,IDLSTA:',LISTA,IDLSTA
380         WRITE (LUPRI,*) 'LISTB,IDLSTB:',LISTB,IDLSTB
381         WRITE (LUPRI,*) 'ISYMA,ISYMB,ISYMAB:',ISYMA,ISYMB,ISYMAB
382         CALL FLSHFO(LUPRI)
383        END IF
384C
385C
386        IF (NYQMMM) THEN
387
388          RSPTYP='K'
389          CALL CCMM_QRTRANSFORMER(WORK(KTHETA1),WORK(KTHETA2),ISYMAB,
390     *                         LISTB,IDLSTB,ISYMB,
391     *                         LISTA,IDLSTA,ISYMA,
392     *                         MODELW,RSPTYP,WORK(KEND1),LWRK1)
393C
394        ELSE IF ((.NOT. NYQMMM) .AND. (.NOT. USE_PELIB())) THEN
395          KTGB  = KEND1
396          KEND2 = KTGB + N2BST(ISYMB)
397          LWRK2 = LWORK   - KEND2
398          IF (LWRK2.LT.0) CALL QUIT( 'Too little work in CC_KMAT, 1')
399C
400          CALL DZERO(WORK(KTGB),N2BST(ISYMB))
401C
402C------------------------------------------------
403C       Trial vector (B left) one excitation part
404C------------------------------------------------
405C
406          KT1AMPB = KEND2
407          KEND3   = KT1AMPB +  NT1AM(ISYMB)
408          LWRK3   = LWORK -  KEND3
409          IF (LWRK3 .LT. 0) THEN
410            CALL QUIT('Insuff. work in CC_KMAT 2')
411          END IF
412          CALL DZERO(WORK(KT1AMPB),NT1AM(ISYMB))
413C
414          IOPT = 1
415          CALL CC_RDRSP(LISTB,IDLSTB,ISYMB,IOPT,MODEL,
416     *                WORK(KT1AMPB),WORK(KDUM))
417C
418          IF (.NOT. CCMM ) CALL CCSL_TGB(WORK(KT1AMPB),ISYMB,
419     *                                 LISTB,IDLSTB,WORK(KTGB),
420     *                                 'XI',MODEL,
421     *                                 WORK(KEND3),LWRK3)
422
423C
424          IF (CCMM) CALL CCMM_TGB(WORK(KT1AMPB),ISYMB,
425     *                          LISTB,IDLSTB,WORK(KTGB),
426     *                         'XI',MODEL,
427     *                          WORK(KEND3),LWRK3)
428
429          NAMPF   = NT1AM(ISYMAB) + NT2AM(ISYMAB)
430C
431          KETA    = KEND2
432          KEND3   = KETA + NAMPF
433          LWRK3   = LWORK   - KEND3
434          IF (LWRK3.LT.0) CALL QUIT( 'Too little work in CC_KMAT, 3')
435          CALL DZERO(WORK(KETA),NAMPF)
436C
437          IF ( LOCDBG ) THEN
438            TAL1= DDOT(N2BST(ISYMAB),WORK(KTGB),1,WORK(KTGB),1)
439            WRITE(LUPRI,*) 'Printing norm2 TGCMO in KMAT: ', TAL1
440          END IF
441
442          LABEL = 'GIVE INT'
443          CALL CC_ETAC(ISYMB,LABEL,WORK(KETA),
444     *                 LISTA,IDLSTA,0,WORK(KTGB),WORK(KEND3),LWRK3)
445
446C
447          KETA1   = KETA
448          KETA2   = KETA + NT1AM(ISYMAB)
449
450          IF (LOCDBG) THEN
451            TAL1= DDOT(NT1AM(ISYMAB),WORK(KETA1),1,WORK(KETA1),1)
452            TAL2= DDOT(NT2AM(ISYMAB),WORK(KETA2),1,
453     *              WORK(KETA2),1)
454            WRITE(LUPRI,*) 'Printing first contribution.
455     &                     Norm2 of singles: ', TAL1,
456     &                    'Norm2 of doubles: ', TAL2
457          END IF
458C
459          LSAME  =  (LISTA.EQ.LISTB  .AND. IDLSTA.EQ.IDLSTB)
460          IF (LSAME) THEN
461            FACTSLV = TWO
462          ELSE
463            FACTSLV = ONE
464          END IF
465C
466          CALL DAXPY(NT1AM(ISYMAB),FACTSLV,WORK(KETA1),1,
467     *               WORK(KTHETA1),1)
468          CALL DAXPY(NT2AM(ISYMAB),FACTSLV,WORK(KETA2),1,
469     *               WORK(KTHETA2),1)
470C
471C------------------------------------------------------------------
472C
473          IF (.NOT. (LSAME)) THEN
474C
475            KTGA  = KEND1
476            KEND2 = KTGA + N2BST(ISYMA)
477            LWRK2 = LWORK   - KEND2
478            IF (LWRK2.LT.0) CALL QUIT( 'Too little work in CC_KMAT, 4')
479            CALL DZERO(WORK(KTGA),N2BST(ISYMA))
480C
481C------------------------------------------------
482C         Trial vector (A left) one excitation part
483C------------------------------------------------
484C
485            KT1AMPA = KEND2
486            KEND3   = KT1AMPA +  NT1AM(ISYMA)
487            LWRK3   = LWORK -  KEND3
488            IF (LWRK3 .LT. 0) THEN
489              CALL QUIT('Insuff. work in CC_KMAT 5')
490            END IF
491            CALL DZERO(WORK(KT1AMPA),NT1AM(ISYMA))
492C
493            IOPT = 1
494            CALL CC_RDRSP(LISTA,IDLSTA,ISYMA,IOPT,MODEL,
495     *                  WORK(KT1AMPA),WORK(KDUM))
496C
497            IF (.NOT. CCMM) CALL CCSL_TGB(WORK(KT1AMPA),ISYMA,LISTA,
498     *                    IDLSTA,WORK(KTGA),'XI',
499     *                    MODEL,WORK(KEND3),LWRK3)
500C
501            IF (CCMM) CALL CCMM_TGB(WORK(KT1AMPA),ISYMA,LISTA,
502     *                   IDLSTA,WORK(KTGA),'XI',
503     *                   MODEL,WORK(KEND3),LWRK3)
504C
505            NAMPF   = NT1AM(ISYMAB) + NT2AM(ISYMAB)
506C
507            KETA    = KEND2
508            KEND3   = KETA + NAMPF
509            LWRK3   = LWORK   - KEND3
510            IF (LWRK3.LT.0) CALL QUIT( 'Too little work in CC_KMAT, 6')
511            CALL DZERO(WORK(KETA),NAMPF)
512C
513            LABEL = 'GIVE INT'
514            CALL CC_ETAC(ISYMA,LABEL,WORK(KETA),
515     *                   LISTB,IDLSTB,0,WORK(KTGA),WORK(KEND3),LWRK3)
516C
517            KETA1   = KETA
518            KETA2   = KETA + NT1AM(ISYMAB)
519
520            IF (LOCDBG) THEN
521              TAL1= DDOT(NT1AM(ISYMAB),WORK(KETA1),1,WORK(KETA1),1)
522              TAL2= DDOT(NT2AM(ISYMAB),WORK(KETA2),1,
523     *              WORK(KETA2),1)
524              WRITE(LUPRI,*) 'Printing second contribution.
525     &                     Norm2 of singles: ', TAL1,
526     &                    'Norm2 of doubles: ', TAL2
527            END IF
528C
529            CALL DAXPY(NT1AM(ISYMAB),ONE,WORK(KETA1),1,WORK(KTHETA1),1)
530            CALL DAXPY(NT2AM(ISYMAB),ONE,WORK(KETA2),1,WORK(KTHETA2),1)
531C
532          END IF
533
534        END IF ! NYQMMM
535C
536        IF (USE_PELIB()) THEN
537          RSPTYP='K'
538          CALL PELIB_IFC_QRTRANSFORMER(WORK(KTHETA1),WORK(KTHETA2),
539     &                   ISYMAB,LISTB,IDLSTB,ISYMB,LISTA,IDLSTA,ISYMA,
540     &                   MODELW,RSPTYP,WORK(KEND1),LWRK1)
541        END IF
542
543*---------------------------------------------------------------------*
544* write result vector to output:
545*---------------------------------------------------------------------*
546
547        IF (IOPTRES .EQ. 0  .OR. IOPTRES .EQ. 1) THEN
548
549*       write to a common direct access file,
550*       store start address in IBTRAN(3,ITRAN)
551
552        IBTRAN(3,ITRAN) = IADRTH
553
554        CALL PUTWA2(LUBMAT,FILBMA,WORK(KTHETA1),IADRTH,NT1AM(ISYMAB))
555        IADRTH = IADRTH + NT1AM(ISYMAB)
556
557        IF (.NOT.CCS) THEN
558          CALL PUTWA2(LUBMAT,FILBMA,WORK(KTHETA2),IADRTH,NT2AM(ISYMAB))
559          IADRTH = IADRTH + NT2AM(ISYMAB)
560        END IF
561
562        IF (LOCDBG) THEN
563         WRITE (LUPRI,*) 'K matrix transformation nb. ',ITRAN,
564     &          ' saved on file.'
565         WRITE (LUPRI,*) 'ADRESS, LENGTH:',
566     &        IBTRAN(3,ITRAN),IADRTH-IBTRAN(3,ITRAN)
567         XNORM = DDOT(NT1AM(ISYMAB),WORK(KTHETA1),1,WORK(KTHETA1),1)
568         IF (.NOT.CCS) XNORM = XNORM +
569     &           DDOT(NT2AM(ISYMAB),WORK(KTHETA2),1,WORK(KTHETA2),1)
570         WRITE (LUPRI,*) 'Norm:', XNORM
571
572         Call AROUND('K matrix transformation written to file:')
573         Call CC_PRP(WORK(KTHETA1),WORK(KTHETA2),ISYMAB,1,1)
574        END IF
575
576      ELSE IF ( IOPTRES .EQ. 3 .OR. IOPTRES .EQ. 4 ) THEN
577
578*        write to a sequential file by a call to CC_WRRSP/CC_WARSP,
579*        use FILBMA as LIST type and IBTRAN(3,ITRAN) as index
580         KTHETA0 = -999999
581         IF (IOPTRES.EQ.3) THEN
582           CALL CC_WRRSP(FILBMA,IBTRAN(3,ITRAN),ISYMAB,IOPTW,MODELW,
583     &                   WORK(KTHETA0),WORK(KTHETA1),WORK(KTHETA2),
584     &                   WORK(KEND1),LWRK1)
585         ELSE IF (IOPTRES.EQ.4) THEN
586           CALL CC_WARSP(FILBMA,IBTRAN(3,ITRAN),ISYMAB,IOPTW,MODELW,
587     &                   WORK(KTHETA0),WORK(KTHETA1),WORK(KTHETA2),
588     &                   WORK(KEND1),LWRK1)
589         END IF
590
591         IF (LOCDBG) THEN
592           WRITE (LUPRI,*) 'Write K * ',LISTA,' * ',LISTB,
593     &              ' transformation',
594     &              ' as ',FILBMA,' type vector to file.'
595           WRITE (LUPRI,*) 'index of inp. A vector:',IBTRAN(1,ITRAN)
596           WRITE (LUPRI,*) 'index of inp. B vector:',IBTRAN(2,ITRAN)
597           WRITE (LUPRI,*) 'index of result vector:',IBTRAN(3,ITRAN)
598           LEN = NT1AM(ISYMAB) + NT2AM(ISYMAB)
599           IF (CCS) LEN = NT1AM(ISYMAB)
600           XNORM = DDOT(LEN,WORK(KTHETA1),1,WORK(KTHETA1),1)
601           WRITE (LUPRI,*) 'norm^2 of result vector:',XNORM
602         END IF
603      ELSE IF (IOPTRES.EQ.5) THEN
604         IF (.NOT.CCS) CALL CCLR_DIASCL(WORK(KTHETA2),TWO,ISYMAB)
605         CALL CCDOTRSP(IBDOTS,BCONS,IOPTW,FILBMA,ITRAN,NBTRAN,MXVEC,
606     &                 WORK(KTHETA1),WORK(KTHETA2),ISYMAB,
607     &                 WORK(KEND1),LWRK1)
608      ELSE
609        CALL QUIT('Illegal value for IOPTRES in CC_KMAT.')
610      END IF
611
612      TIMTRN = SECOND() - TIMTRN
613
614      IF (IPRINT.GT.0) THEN
615
616         IF (IOPTRES.EQ.5) THEN
617            IVEC = 1
618            DO WHILE (IBDOTS(IVEC,ITRAN).NE.0 .AND. IVEC.LE.MXVEC)
619               IVEC = IVEC + 1
620            END DO
621            WRITE (LUPRI,'(1X,2(A,I5),A,I6,A,F10.2,A)')'| ',IDLSTA,
622     &        '    | ',IDLSTB,'    | ',IVEC-1,'       | ',TIMTRN,'  |'
623         ELSE
624            WRITE (LUPRI,'(1X,2(A,I5),A,I6,A,F10.2,A)') '| ',IDLSTA,
625     &           '    | ',
626     &        IDLSTB,'    | ',IBTRAN(3,ITRAN),'       | ',TIMTRN,'  |'
627         END IF
628
629      END IF
630
631*---------------------------------------------------------------------*
632* End of loop over K matrix transformations
633*---------------------------------------------------------------------*
634      END DO
635      WRITE (LUPRI,'(1X,A1,50("="),A1,//)') '+','+'
636
637*---------------------------------------------------------------------*
638* if IOPTRES=1 and enough work space available, read result
639* vectors back into memory:
640*---------------------------------------------------------------------*
641
642* check size of work space:
643      IF (IOPTRES .EQ. 1) THEN
644        LENALL = IADRTH-1
645        IF (LENALL .GT. LWORK) IOPTRES = 0
646      END IF
647
648* read the result vectors back into memory:
649      IF (IOPTRES .EQ. 1) THEN
650
651        CALL GETWA2(LUBMAT,FILBMA,WORK(1),1,LENALL)
652
653        IF (LOCDBG) THEN
654          DO ITRAN = 1, NBTRAN
655            IF (ITRAN.LT.NBTRAN) THEN
656              LEN     = IBTRAN(3,ITRAN+1)-IBTRAN(3,ITRAN)
657            ELSE
658              LEN     = IADRTH-IBTRAN(3,NBTRAN)
659            END IF
660            KTHETA1 = IBTRAN(3,ITRAN)
661            XNORM   = DDOT(LEN, WORK(KTHETA1),1, WORK(KTHETA1),1)
662            WRITE (LUPRI,*) 'Read K matrix transformation nb. ',NBTRAN
663            WRITE (LUPRI,*) 'Adress, length, NORM:',IBTRAN(3,NBTRAN),
664     &                      LEN,XNORM
665          END DO
666          CALL FLSHFO(LUPRI)
667        END IF
668      END IF
669
670*---------------------------------------------------------------------*
671* close K matrix file, print timings & return
672*---------------------------------------------------------------------*
673
674      IF (IOPTRES.EQ.0 ) THEN
675        CALL WCLOSE2(LUBMAT, FILBMA, 'KEEP')
676      ELSE IF (IOPTRES.EQ.1) THEN
677        CALL WCLOSE2(LUBMAT, FILBMA, 'DELETE')
678      ELSE IF (IOPTRES.EQ.3 .OR. IOPTRES.EQ.4 .OR. IOPTRES.EQ.5) THEN
679        CONTINUE
680      ELSE
681        CALL QUIT('Illegal value of IOPTRES in CC_KMAT.')
682      END IF
683
684*=====================================================================*
685
686      RETURN
687      END
688*=====================================================================*
689*            END OF SUBROUTINE CC_KMAT
690*=====================================================================*
691
692