1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
|
C Copyright 1981-2007 ECMWF
C
C Licensed under the GNU Lesser General Public License which
C incorporates the terms and conditions of version 3 of the GNU
C General Public License.
C See LICENSE and gpl-3.0.txt for details.
C
INTEGER FUNCTION IGTRAN (PMAT, KINCOL, KINROW, KMOVE, KLEN,
1 KPR, KERR)
C
C---->
C**** *IGTRAN*
C
C PURPOSE
C _______
C
C This routine transposes a rectangular matrix.
C
C INTERFACE
C _________
C
C IERR = IGTRAN (PMAT, KINCOL, KINROW, KMOVE, KLEN, KPR, KERR)
C
C Input parameters
C ________________
C
C PMAT - The input matrix of length KINCOL * KINROW.
C
C KINCOL - The length of the first dimension on entry and
C second dimension on exit.
C
C KINROW - The length of the second dimension on entry and
C first dimension on exit.
C
C KMOVE - Work array used to store information to speed up
C the process.
C
C KLEN - Length of array KMOVE. The recommended length is
C (KINCOL + KINROW) / 2.
C
C KPR - The debug print switch.
C 0 , No debugging output.
C 1 , Produce debugging output.
C
C KERR - The error control flag.
C (No longer used, kept for backward compatibility.
C Used to be used as follows:
C -ve, No error message. Return error code.
C 0 , Hard failure with error message.
C +ve, Print error message. Return error code.)
C
C Output parameters
C ________________
C
C PMAT - The transposed matrix.
C
C Return value
C ____________
C
C The error indicator (INTEGER).
C
C 7401 KLEN was less than 1.
C 7402 A failure during transposition (should never happen).
C
C
C Common block usage
C __________________
C
C None
C
C EXTERNALS
C _________
C
C INTLOG - Generate log messages.
C
C
C METHOD
C ______
C
C This algorithm uses the cyclic structure of transposition to
C perform a transposition in place with a minimum amount of work
C storage.
C
C REFERENCE
C _________
C
C Esko G. Cate and David W. Twigg Analysis of In-Situ
C Transposition
C CACM Algorithm 513
C
C
C COMMENTS
C ________
C
C None.
C
C
C AUTHOR
C ______
C
C K. Fielding *ECMWF* Jan 1994
C
C
C MODIFICATIONS
C _____________
C
C J.D.Chambers ECMWF Sept 1995
C
C----<
C _______________________________________________________
C
IMPLICIT NONE
C
#include "parim.h"
C
C Function arguments
C
INTEGER KINCOL, KINROW, KLEN, KPR, KERR
INTEGER KMOVE(KLEN)
REAL PMAT (KINCOL * KINROW)
C
C Local variables
INTEGER ICOUNT, IRMULC, IR0, IR1, IR2, IPOINT, IKMI, IP1, IP2,
X IRMCM1, IP1C, IP2C, IROWP1, ICOLM1, ISTART
INTEGER JST, JINROW, JINCOL
REAL ZMIP1, ZMIP1C, ZTEMP
C
C Parameters
INTEGER JPROUTINE
PARAMETER (JPROUTINE = 7400)
C
C _______________________________________________________
C
C* Section 1. Initialisation
C _______________________________________________________
C
100 CONTINUE
C
IF (KPR .GE. 1) CALL INTLOG(JP_DEBUG,'IGTRAN: Section 1.',JPQUIET)
C
IGTRAN = 0
C
IF (KPR .GE. 1) THEN
CALL INTLOG(JP_DEBUG,'IGTRAN: Input parameters.',JPQUIET)
CALL INTLOG(JP_DEBUG,'IGTRAN: 1st matrix dimension = ',KINCOL)
CALL INTLOG(JP_DEBUG,'IGTRAN: 2nd matrix dimension = ',KINROW)
CALL INTLOG(JP_DEBUG,'IGTRAN: Work array length = ',KLEN)
CALL INTLOG(JP_DEBUG,
X 'IGTRAN: Recommended length = ',(KINCOL + KINROW) / 2)
ENDIF
C
IF (KINCOL .LE. 1 .OR. KINROW .LE. 1) GO TO 900
C
C _______________________________________________________
C
C* Section 2. Rectangular transposition setup
C _______________________________________________________
C
200 CONTINUE
C
IF (KPR .GE. 1) CALL INTLOG(JP_DEBUG,'IGTRAN: Section 2.',JPQUIET)
C
IF (KINCOL .NE. KINROW) THEN
C
IF (KLEN .LT. 1) THEN
IGTRAN = JPROUTINE + 1
CALL INTLOG(JP_ERROR,'IGTRAN: Work array size = ',KLEN)
CALL INTLOG(JP_ERROR,'IGTRAN: Must be at least 1.',JPQUIET)
GO TO 900
ENDIF
C
ICOUNT = 2
IRMULC = KINCOL * KINROW
IRMCM1 = IRMULC - 1
C
DO 210 JST = 1, KLEN
KMOVE (JST) = 0
210 CONTINUE
C
IF (KINCOL .GT. 2 .AND. KINROW .GT. 2) THEN
C
C Calculate the number of fixed points, Euclids algorithm
C for GCD (m - 1, n - 1)
C
IR2 = KINCOL - 1
IR1 = KINROW - 1
C
220 CONTINUE
IR0 = MOD (IR2, IR1)
IR2 = IR1
IR1 = IR0
IF (IR0 .NE. 0) GO TO 220
C
ICOUNT = ICOUNT + IR2 - 1
C
ENDIF
C
C Set initial values for search
C
ISTART = 1
IPOINT = KINCOL
C
C At least one loop must be rearranged so branch into loop
C
GO TO 330
C
C _______________________________________________________
C
C* Section 3. Rectangular transposition main loop
C _______________________________________________________
C
310 CONTINUE
C
C Search for loops to rearrange
C
IKMI = IRMCM1 - ISTART
ISTART = ISTART + 1
C
IF (ISTART .GT. IKMI) THEN
IGTRAN = JPROUTINE + 2
CALL INTLOG(JP_ERROR,
X 'IGTRAN: Fail during transposition.',JPQUIET)
GO TO 900
ENDIF
C
IPOINT = IPOINT + KINCOL
C
IF (IPOINT .GT. IRMCM1) IPOINT = IPOINT - IRMCM1
C
IP2 = IPOINT
C
IF (ISTART .EQ. IP2) GO TO 310
C
IF (ISTART .GT. KLEN) THEN
C
320 CONTINUE
C
C Loop exit condition
C
IF (IP2 .LE. ISTART .OR. IP2 .GE. IKMI) THEN
IF (IP2 .NE. ISTART) THEN
GO TO 310
ELSE
GO TO 330
ENDIF
ENDIF
IP1 = IP2
IP2 = KINCOL * IP1 - IRMCM1 * (IP1 / KINROW)
GO TO 320
ENDIF
C
IF (KMOVE (ISTART) .NE. 0) GO TO 310
C
C Rearrange the elements of a loop and its companion loop
C
C Entry into loop on first pass
C
330 CONTINUE
C
IP1 = ISTART
IKMI = IRMCM1 - ISTART
ZMIP1 = PMAT (IP1 + 1)
IP1C = IKMI
ZMIP1C = PMAT (IP1C + 1)
C
340 CONTINUE
IP2 = KINCOL * IP1 - IRMCM1 * (IP1 / KINROW)
IP2C = IRMCM1 - IP2
C
IF (IP1 .LE. KLEN) KMOVE (IP1) = 2
IF (IP1C .LE. KLEN) KMOVE (IP1C) = 2
C
ICOUNT = ICOUNT + 2
C
C Loop exit conditions
C
IF (IP2 .EQ. ISTART) GO TO 360
IF (IP2 .EQ. IKMI) GO TO 350
C
PMAT (IP1 + 1) = PMAT(IP2 + 1)
PMAT (IP1C + 1) = PMAT(IP2C + 1)
IP1 = IP2
IP1C = IP2C
C
GO TO 340
C
350 CONTINUE
C
ZTEMP = ZMIP1
ZMIP1 = ZMIP1C
ZMIP1C = ZTEMP
C
360 CONTINUE
C
PMAT (IP1 + 1) = ZMIP1
PMAT (IP1C + 1) = ZMIP1C
C
IF (ICOUNT .LT. IRMULC) GO TO 310
C
ELSE
C
C _______________________________________________________
C
C* Section 4. Square transposition
C _______________________________________________________
C
400 CONTINUE
C
IF (KPR.GE.1) CALL INTLOG(JP_DEBUG,'IGTRAN: Section 4.',JPQUIET)
C
C Square matrix so exchange elements a(i,j) and a(j,i)
C
ICOLM1 = KINCOL - 1
C
DO 420 JINCOL = 1, ICOLM1
C
IROWP1 = JINCOL + 1
C
DO 410 JINROW = IROWP1, KINROW
C
IP1 = JINCOL + (JINROW - 1) * KINCOL
IP2 = JINROW + (JINCOL - 1) * KINROW
C
ZTEMP = PMAT (IP1)
PMAT (IP1) = PMAT (IP2)
PMAT (IP2) = ZTEMP
C
410 CONTINUE
420 CONTINUE
ENDIF
C
C _______________________________________________________
C
C* Section 9. Return to calling routine. Format statements
C _______________________________________________________
C
900 CONTINUE
C
IF (KPR .GE. 1) CALL INTLOG(JP_DEBUG,'IGTRAN: Section 9.',JPQUIET)
C
RETURN
END
|
