Line data Source code
1 : /* pack_gp.f -- translated by f2c (version 20031025).
2 : You must link the resulting object file with libf2c:
3 : on Microsoft Windows system, link with libf2c.lib;
4 : on Linux or Unix systems, link with .../path/to/libf2c.a -lm
5 : or, if you install libf2c.a in a standard place, with -lf2c -lm
6 : -- in that order, at the end of the command line, as in
7 : cc *.o -lf2c -lm
8 : Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
9 :
10 : http://www.netlib.org/f2c/libf2c.zip
11 : */
12 :
13 : /*#include "f2c.h"*/
14 : #include "cpl_port.h"
15 : #include <limits.h>
16 : #include <stdlib.h>
17 : #include "grib2.h"
18 : typedef g2int logical;
19 : #define TRUE_ (1)
20 : #define FALSE_ (0)
21 :
22 13 : /* Subroutine */ int pack_gp(integer *kfildo, integer *ic, integer *nxy,
23 : integer *is523, integer *minpk, integer *inc, integer *missp, integer
24 : *misss, integer *jmin, integer *jmax, integer *lbit, integer *nov,
25 : integer *ndg, integer *lx, integer *ibit, integer *jbit, integer *
26 : kbit, integer *novref, integer *lbitref, integer *ier)
27 : {
28 : /* Initialized data */
29 :
30 13 : const integer mallow = 1073741825; /* MALLOW=2**30+1 */
31 13 : integer ifeed = 12;
32 13 : integer ifirst = 0;
33 :
34 : /* System generated locals */
35 : integer i__1, i__2, i__3;
36 :
37 : /* Local variables */
38 : integer j, k, l;
39 : logical adda;
40 13 : integer ired, kinc, mina, maxa, minb = 0, maxb = 0, minc = 0, maxc = 0, ibxx2[31];
41 : char cfeed[1];
42 13 : integer nenda, nendb = 0, ibita, ibitb = 0, minak, minbk = 0, maxak, maxbk = 0,
43 13 : minck = 0, maxck, nouta = 0, lmiss, itest = 0, nount = 0;
44 : extern /* Subroutine */ int reduce(integer *, integer *, integer *,
45 : integer *, integer *, integer *, integer *, integer *, integer *,
46 : integer *, integer *, integer *, integer *);
47 13 : integer ibitbs = 0, mislla, misllb = 0, misllc = 0, iersav = 0, lminpk, ktotal,
48 13 : kounta, kountb = 0, kstart, mstart = 0, mintst = 0, maxtst = 0,
49 13 : kounts = 0, mintstk = 0, maxtstk = 0;
50 : integer *misslx;
51 :
52 :
53 : /* FEBRUARY 1994 GLAHN TDL MOS-2000 */
54 : /* JUNE 1995 GLAHN MODIFIED FOR LMISS ERROR. */
55 : /* JULY 1996 GLAHN ADDED MISSS */
56 : /* FEBRUARY 1997 GLAHN REMOVED 4 REDUNDANT TESTS FOR */
57 : /* MISSP.EQ.0; INSERTED A TEST TO BETTER */
58 : /* HANDLE A STRING OF 9999'S */
59 : /* FEBRUARY 1997 GLAHN ADDED LOOPS TO ELIMINATE TEST FOR */
60 : /* MISSS WHEN MISSS = 0 */
61 : /* MARCH 1997 GLAHN CORRECTED FOR SECONDARY MISSING VALUE */
62 : /* MARCH 1997 GLAHN CORRECTED FOR USE OF LOCAL VALUE */
63 : /* OF MINPK */
64 : /* MARCH 1997 GLAHN CORRECTED FOR SECONDARY MISSING VALUE */
65 : /* MARCH 1997 GLAHN CHANGED CALCULATING NUMBER OF BITS */
66 : /* THROUGH EXPONENTS TO AN ARRAY (IMPROVED */
67 : /* OVERALL PACKING PERFORMANCE BY ABOUT */
68 : /* 35 PERCENT!). ALLOWED 0 BITS FOR */
69 : /* PACKING JMIN( ), LBIT( ), AND NOV( ). */
70 : /* MAY 1997 GLAHN A NUMBER OF CHANGES FOR EFFICIENCY. */
71 : /* MOD FUNCTIONS ELIMINATED AND ONE */
72 : /* IFTHEN ADDED. JOUNT REMOVED. */
73 : /* RECOMPUTATION OF BITS NOT MADE UNLESS */
74 : /* NECESSARY AFTER MOVING POINTS FROM */
75 : /* ONE GROUP TO ANOTHER. NENDB ADJUSTED */
76 : /* TO ELIMINATE POSSIBILITY OF VERY */
77 : /* SMALL GROUP AT THE END. */
78 : /* ABOUT 8 PERCENT IMPROVEMENT IN */
79 : /* OVERALL PACKING. ISKIPA REMOVED; */
80 : /* THERE IS ALWAYS A GROUP B THAT CAN */
81 : /* BECOME GROUP A. CONTROL ON SIZE */
82 : /* OF GROUP B (STATEMENT BELOW 150) */
83 : /* ADDED. ADDED ADDA, AND USE */
84 : /* OF GE AND LE INSTEAD OF GT AND LT */
85 : /* IN LOOPS BETWEEN 150 AND 160. */
86 : /* IBITBS ADDED TO SHORTEN TRIPS */
87 : /* THROUGH LOOP. */
88 : /* MARCH 2000 GLAHN MODIFIED FOR GRIB2; CHANGED NAME FROM */
89 : /* PACKGP */
90 : /* JANUARY 2001 GLAHN COMMENTS; IER = 706 SUBSTITUTED FOR */
91 : /* STOPS; ADDED RETURN1; REMOVED STATEMENT */
92 : /* NUMBER 110; ADDED IER AND * RETURN */
93 : /* NOVEMBER 2001 GLAHN CHANGED SOME DIAGNOSTIC FORMATS TO */
94 : /* ALLOW PRINTING LARGER NUMBERS */
95 : /* NOVEMBER 2001 GLAHN ADDED MISSLX( ) TO PUT MAXIMUM VALUE */
96 : /* INTO JMIN( ) WHEN ALL VALUES MISSING */
97 : /* TO AGREE WITH GRIB STANDARD. */
98 : /* NOVEMBER 2001 GLAHN CHANGED TWO TESTS ON MISSP AND MISSS */
99 : /* EQ 0 TO TESTS ON IS523. HOWEVER, */
100 : /* MISSP AND MISSS CANNOT IN GENERAL BE */
101 : /* = 0. */
102 : /* NOVEMBER 2001 GLAHN ADDED CALL TO REDUCE; DEFINED ITEST */
103 : /* BEFORE LOOPS TO REDUCE COMPUTATION; */
104 : /* STARTED LARGE GROUP WHEN ALL SAME */
105 : /* VALUE */
106 : /* DECEMBER 2001 GLAHN MODIFIED AND ADDED A FEW COMMENTS */
107 : /* JANUARY 2002 GLAHN REMOVED LOOP BEFORE 150 TO DETERMINE */
108 : /* A GROUP OF ALL SAME VALUE */
109 : /* JANUARY 2002 GLAHN CHANGED MALLOW FROM 9999999 TO 2**30+1, */
110 : /* AND MADE IT A PARAMETER */
111 : /* MARCH 2002 GLAHN ADDED NON FATAL IER = 716, 717; */
112 : /* REMOVED NENDB=NXY ABOVE 150; */
113 : /* ADDED IERSAV=0; COMMENTS */
114 :
115 : /* PURPOSE */
116 : /* DETERMINES GROUPS OF VARIABLE SIZE, BUT AT LEAST OF */
117 : /* SIZE MINPK, THE ASSOCIATED MAX (JMAX( )) AND MIN (JMIN( )), */
118 : /* THE NUMBER OF BITS NECESSARY TO HOLD THE VALUES IN EACH */
119 : /* GROUP (LBIT( )), THE NUMBER OF VALUES IN EACH GROUP */
120 : /* (NOV( )), THE NUMBER OF BITS NECESSARY TO PACK THE JMIN( ) */
121 : /* VALUES (IBIT), THE NUMBER OF BITS NECESSARY TO PACK THE */
122 : /* LBIT( ) VALUES (JBIT), AND THE NUMBER OF BITS NECESSARY */
123 : /* TO PACK THE NOV( ) VALUES (KBIT). THE ROUTINE IS DESIGNED */
124 : /* TO DETERMINE THE GROUPS SUCH THAT A SMALL NUMBER OF BITS */
125 : /* IS NECESSARY TO PACK THE DATA WITHOUT EXCESSIVE */
126 : /* COMPUTATIONS. IF ALL VALUES IN THE GROUP ARE ZERO, THE */
127 : /* NUMBER OF BITS TO USE IN PACKING IS DEFINED AS ZERO WHEN */
128 : /* THERE CAN BE NO MISSING VALUES; WHEN THERE CAN BE MISSING */
129 : /* VALUES, THE NUMBER OF BITS MUST BE AT LEAST 1 TO HAVE */
130 : /* THE CAPABILITY TO RECOGNIZE THE MISSING VALUE. HOWEVER, */
131 : /* IF ALL VALUES IN A GROUP ARE MISSING, THE NUMBER OF BITS */
132 : /* NEEDED IS 0, AND THE UNPACKER RECOGNIZES THIS. */
133 : /* ALL VARIABLES ARE INTEGER. EVEN THOUGH THE GROUPS ARE */
134 : /* INITIALLY OF SIZE MINPK OR LARGER, AN ADJUSTMENT BETWEEN */
135 : /* TWO GROUPS (THE LOOKBACK PROCEDURE) MAY MAKE A GROUP */
136 : /* SMALLER THAN MINPK. THE CONTROL ON GROUP SIZE IS THAT */
137 : /* THE SUM OF THE SIZES OF THE TWO CONSECUTIVE GROUPS, EACH OF */
138 : /* SIZE MINPK OR LARGER, IS NOT DECREASED. WHEN DETERMINING */
139 : /* THE NUMBER OF BITS NECESSARY FOR PACKING, THE LARGEST */
140 : /* VALUE THAT CAN BE ACCOMMODATED IN, SAY, MBITS, IS */
141 : /* 2**MBITS-1; THIS LARGEST VALUE (AND THE NEXT SMALLEST */
142 : /* VALUE) IS RESERVED FOR THE MISSING VALUE INDICATOR (ONLY) */
143 : /* WHEN IS523 NE 0. IF THE DIMENSION NDG */
144 : /* IS NOT LARGE ENOUGH TO HOLD ALL THE GROUPS, THE LOCAL VALUE */
145 : /* OF MINPK IS INCREASED BY 50 PERCENT. THIS IS REPEATED */
146 : /* UNTIL NDG WILL SUFFICE. A DIAGNOSTIC IS PRINTED WHENEVER */
147 : /* THIS HAPPENS, WHICH SHOULD BE VERY RARELY. IF IT HAPPENS */
148 : /* OFTEN, NDG IN SUBROUTINE PACK SHOULD BE INCREASED AND */
149 : /* A CORRESPONDING INCREASE IN SUBROUTINE UNPACK MADE. */
150 : /* CONSIDERABLE CODE IS PROVIDED SO THAT NO MORE CHECKING */
151 : /* FOR MISSING VALUES WITHIN LOOPS IS DONE THAN NECESSARY; */
152 : /* THE ADDED EFFICIENCY OF THIS IS RELATIVELY MINOR, */
153 : /* BUT DOES NO HARM. FOR GRIB2, THE REFERENCE VALUE FOR */
154 : /* THE LENGTH OF GROUPS IN NOV( ) AND FOR THE NUMBER OF */
155 : /* BITS NECESSARY TO PACK GROUP VALUES ARE DETERMINED, */
156 : /* AND SUBTRACTED BEFORE JBIT AND KBIT ARE DETERMINED. */
157 :
158 : /* WHEN 1 OR MORE GROUPS ARE LARGE COMPARED TO THE OTHERS, */
159 : /* THE WIDTH OF ALL GROUPS MUST BE AS LARGE AS THE LARGEST. */
160 : /* A SUBROUTINE REDUCE BREAKS UP LARGE GROUPS INTO 2 OR */
161 : /* MORE TO REDUCE TOTAL BITS REQUIRED. IF REDUCE SHOULD */
162 : /* ABORT, PACK_GP WILL BE EXECUTED AGAIN WITHOUT THE CALL */
163 : /* TO REDUCE. */
164 :
165 : /* DATA SET USE */
166 : /* KFILDO - UNIT NUMBER FOR OUTPUT (PRINT) FILE. (OUTPUT) */
167 :
168 : /* VARIABLES IN CALL SEQUENCE */
169 : /* KFILDO = UNIT NUMBER FOR OUTPUT (PRINT) FILE. (INPUT) */
170 : /* IC( ) = ARRAY TO HOLD DATA FOR PACKING. THE VALUES */
171 : /* DO NOT HAVE TO BE POSITIVE AT THIS POINT, BUT */
172 : /* MUST BE IN THE RANGE -2**30 TO +2**30 (THE */
173 : /* THE VALUE OF MALLOW). THESE INTEGER VALUES */
174 : /* WILL BE RETAINED EXACTLY THROUGH PACKING AND */
175 : /* UNPACKING. (INPUT) */
176 : /* NXY = NUMBER OF VALUES IN IC( ). ALSO TREATED */
177 : /* AS ITS DIMENSION. (INPUT) */
178 : /* IS523 = missing value management */
179 : /* 0=data contains no missing values */
180 : /* 1=data contains Primary missing values */
181 : /* 2=data contains Primary and secondary missing values */
182 : /* (INPUT) */
183 : /* MINPK = THE MINIMUM SIZE OF EACH GROUP, EXCEPT POSSIBLY */
184 : /* THE LAST ONE. (INPUT) */
185 : /* INC = THE NUMBER OF VALUES TO ADD TO AN ALREADY */
186 : /* EXISTING GROUP IN DETERMINING WHETHER OR NOT */
187 : /* TO START A NEW GROUP. IDEALLY, THIS WOULD BE */
188 : /* 1, BUT EACH TIME INC VALUES ARE ATTEMPTED, THE */
189 : /* MAX AND MIN OF THE NEXT MINPK VALUES MUST BE */
190 : /* FOUND. THIS IS "A LOOP WITHIN A LOOP," AND */
191 : /* A SLIGHTLY LARGER VALUE MAY GIVE ABOUT AS GOOD */
192 : /* RESULTS WITH SLIGHTLY LESS COMPUTATIONAL TIME. */
193 : /* IF INC IS LE 0, 1 IS USED, AND A DIAGNOSTIC IS */
194 : /* OUTPUT. NOTE: IT IS EXPECTED THAT INC WILL */
195 : /* EQUAL 1. THE CODE USES INC PRIMARILY IN THE */
196 : /* LOOPS STARTING AT STATEMENT 180. IF INC */
197 : /* WERE 1, THERE WOULD NOT NEED TO BE LOOPS */
198 : /* AS SUCH. HOWEVER, KINC (THE LOCAL VALUE OF */
199 : /* INC) IS SET GE 1 WHEN NEAR THE END OF THE DATA */
200 : /* TO FORESTALL A VERY SMALL GROUP AT THE END. */
201 : /* (INPUT) */
202 : /* MISSP = WHEN MISSING POINTS CAN BE PRESENT IN THE DATA, */
203 : /* THEY WILL HAVE THE VALUE MISSP OR MISSS. */
204 : /* MISSP IS THE PRIMARY MISSING VALUE AND MISSS */
205 : /* IS THE SECONDARY MISSING VALUE . THESE MUST */
206 : /* NOT BE VALUES THAT WOULD OCCUR WITH SUBTRACTING */
207 : /* THE MINIMUM (REFERENCE) VALUE OR SCALING. */
208 : /* FOR EXAMPLE, MISSP = 0 WOULD NOT BE ADVISABLE. */
209 : /* (INPUT) */
210 : /* MISSS = SECONDARY MISSING VALUE INDICATOR (SEE MISSP). */
211 : /* (INPUT) */
212 : /* JMIN(J) = THE MINIMUM OF EACH GROUP (J=1,LX). (OUTPUT) */
213 : /* JMAX(J) = THE MAXIMUM OF EACH GROUP (J=1,LX). THIS IS */
214 : /* NOT REALLY NEEDED, BUT SINCE THE MAX OF EACH */
215 : /* GROUP MUST BE FOUND, SAVING IT HERE IS CHEAP */
216 : /* IN CASE THE USER WANTS IT. (OUTPUT) */
217 : /* LBIT(J) = THE NUMBER OF BITS NECESSARY TO PACK EACH GROUP */
218 : /* (J=1,LX). IT IS ASSUMED THE MINIMUM OF EACH */
219 : /* GROUP WILL BE REMOVED BEFORE PACKING, AND THE */
220 : /* VALUES TO PACK WILL, THEREFORE, ALL BE POSITIVE. */
221 : /* HOWEVER, IC( ) DOES NOT NECESSARILY CONTAIN */
222 : /* ALL POSITIVE VALUES. IF THE OVERALL MINIMUM */
223 : /* HAS BEEN REMOVED (THE USUAL CASE), THEN IC( ) */
224 : /* WILL CONTAIN ONLY POSITIVE VALUES. (OUTPUT) */
225 : /* NOV(J) = THE NUMBER OF VALUES IN EACH GROUP (J=1,LX). */
226 : /* (OUTPUT) */
227 : /* NDG = THE DIMENSION OF JMIN( ), JMAX( ), LBIT( ), AND */
228 : /* NOV( ). (INPUT) */
229 : /* LX = THE NUMBER OF GROUPS DETERMINED. (OUTPUT) */
230 : /* IBIT = THE NUMBER OF BITS NECESSARY TO PACK THE JMIN(J) */
231 : /* VALUES, J=1,LX. (OUTPUT) */
232 : /* JBIT = THE NUMBER OF BITS NECESSARY TO PACK THE LBIT(J) */
233 : /* VALUES, J=1,LX. (OUTPUT) */
234 : /* KBIT = THE NUMBER OF BITS NECESSARY TO PACK THE NOV(J) */
235 : /* VALUES, J=1,LX. (OUTPUT) */
236 : /* NOVREF = REFERENCE VALUE FOR NOV( ). (OUTPUT) */
237 : /* LBITREF = REFERENCE VALUE FOR LBIT( ). (OUTPUT) */
238 : /* IER = ERROR RETURN. */
239 : /* 706 = VALUE WILL NOT PACK IN 30 BITS--FATAL */
240 : /* 714 = ERROR IN REDUCE--NON-FATAL */
241 : /* 715 = NGP NOT LARGE ENOUGH IN REDUCE--NON-FATAL */
242 : /* 716 = MINPK INCEASED--NON-FATAL */
243 : /* 717 = INC SET = 1--NON-FATAL */
244 : /* (OUTPUT) */
245 : /* * = ALTERNATE RETURN WHEN IER NE 0 AND FATAL ERROR. */
246 :
247 : /* INTERNAL VARIABLES */
248 : /* CFEED = CONTAINS THE CHARACTER REPRESENTATION */
249 : /* OF A PRINTER FORM FEED. */
250 : /* IFEED = CONTAINS THE INTEGER VALUE OF A PRINTER */
251 : /* FORM FEED. */
252 : /* KINC = WORKING COPY OF INC. MAY BE MODIFIED. */
253 : /* MINA = MINIMUM VALUE IN GROUP A. */
254 : /* MAXA = MAXIMUM VALUE IN GROUP A. */
255 : /* NENDA = THE PLACE IN IC( ) WHERE GROUP A ENDS. */
256 : /* KSTART = THE PLACE IN IC( ) WHERE GROUP A STARTS. */
257 : /* IBITA = NUMBER OF BITS NEEDED TO HOLD VALUES IN GROUP A. */
258 : /* MINB = MINIMUM VALUE IN GROUP B. */
259 : /* MAXB = MAXIMUM VALUE IN GROUP B. */
260 : /* NENDB = THE PLACE IN IC( ) WHERE GROUP B ENDS. */
261 : /* IBITB = NUMBER OF BITS NEEDED TO HOLD VALUES IN GROUP B. */
262 : /* MINC = MINIMUM VALUE IN GROUP C. */
263 : /* MAXC = MAXIMUM VALUE IN GROUP C. */
264 : /* KTOTAL = COUNT OF NUMBER OF VALUES IN IC( ) PROCESSED. */
265 : /* NOUNT = NUMBER OF VALUES ADDED TO GROUP A. */
266 : /* LMISS = 0 WHEN IS523 = 0. WHEN PACKING INTO A */
267 : /* SPECIFIC NUMBER OF BITS, SAY MBITS, */
268 : /* THE MAXIMUM VALUE THAT CAN BE HANDLED IS */
269 : /* 2**MBITS-1. WHEN IS523 = 1, INDICATING */
270 : /* PRIMARY MISSING VALUES, THIS MAXIMUM VALUE */
271 : /* IS RESERVED TO HOLD THE PRIMARY MISSING VALUE */
272 : /* INDICATOR AND LMISS = 1. WHEN IS523 = 2, */
273 : /* THE VALUE JUST BELOW THE MAXIMUM (I.E., */
274 : /* 2**MBITS-2) IS RESERVED TO HOLD THE SECONDARY */
275 : /* MISSING VALUE INDICATOR AND LMISS = 2. */
276 : /* LMINPK = LOCAL VALUE OF MINPK. THIS WILL BE ADJUSTED */
277 : /* UPWARD WHENEVER NDG IS NOT LARGE ENOUGH TO HOLD */
278 : /* ALL THE GROUPS. */
279 : /* MALLOW = THE LARGEST ALLOWABLE VALUE FOR PACKING. */
280 : /* MISLLA = SET TO 1 WHEN ALL VALUES IN GROUP A ARE MISSING. */
281 : /* THIS IS USED TO DISTINGUISH BETWEEN A REAL */
282 : /* MINIMUM WHEN ALL VALUES ARE NOT MISSING */
283 : /* AND A MINIMUM THAT HAS BEEN SET TO ZERO WHEN */
284 : /* ALL VALUES ARE MISSING. 0 OTHERWISE. */
285 : /* NOTE THAT THIS DOES NOT DISTINGUISH BETWEEN */
286 : /* PRIMARY AND SECONDARY MISSING WHEN SECONDARY */
287 : /* MISSING ARE PRESENT. THIS MEANS THAT */
288 : /* LBIT( ) WILL NOT BE ZERO WITH THE RESULTING */
289 : /* COMPRESSION EFFICIENCY WHEN SECONDARY MISSING */
290 : /* ARE PRESENT. ALSO NOTE THAT A CHECK HAS BEEN */
291 : /* MADE EARLIER TO DETERMINE THAT SECONDARY */
292 : /* MISSING ARE REALLY THERE. */
293 : /* MISLLB = SET TO 1 WHEN ALL VALUES IN GROUP B ARE MISSING. */
294 : /* THIS IS USED TO DISTINGUISH BETWEEN A REAL */
295 : /* MINIMUM WHEN ALL VALUES ARE NOT MISSING */
296 : /* AND A MINIMUM THAT HAS BEEN SET TO ZERO WHEN */
297 : /* ALL VALUES ARE MISSING. 0 OTHERWISE. */
298 : /* MISLLC = PERFORMS THE SAME FUNCTION FOR GROUP C THAT */
299 : /* MISLLA AND MISLLB DO FOR GROUPS B AND C, */
300 : /* RESPECTIVELY. */
301 : /* IBXX2(J) = AN ARRAY THAT WHEN THIS ROUTINE IS FIRST ENTERED */
302 : /* IS SET TO 2**J, J=0,30. IBXX2(30) = 2**30, WHICH */
303 : /* IS THE LARGEST VALUE PACKABLE, BECAUSE 2**31 */
304 : /* IS LARGER THAN THE INTEGER WORD SIZE. */
305 : /* IFIRST = SET BY DATA STATEMENT TO 0. CHANGED TO 1 ON */
306 : /* FIRST */
307 : /* ENTRY WHEN IBXX2( ) IS FILLED. */
308 : /* MINAK = KEEPS TRACK OF THE LOCATION IN IC( ) WHERE THE */
309 : /* MINIMUM VALUE IN GROUP A IS LOCATED. */
310 : /* MAXAK = DOES THE SAME AS MINAK, EXCEPT FOR THE MAXIMUM. */
311 : /* MINBK = THE SAME AS MINAK FOR GROUP B. */
312 : /* MAXBK = THE SAME AS MAXAK FOR GROUP B. */
313 : /* MINCK = THE SAME AS MINAK FOR GROUP C. */
314 : /* MAXCK = THE SAME AS MAXAK FOR GROUP C. */
315 : /* ADDA = KEEPS TRACK WHETHER OR NOT AN ATTEMPT TO ADD */
316 : /* POINTS TO GROUP A WAS MADE. IF SO, THEN ADDA */
317 : /* KEEPS FROM TRYING TO PUT ONE BACK INTO B. */
318 : /* (LOGICAL) */
319 : /* IBITBS = KEEPS CURRENT VALUE IF IBITB SO THAT LOOP */
320 : /* ENDING AT 166 DOESN'T HAVE TO START AT */
321 : /* IBITB = 0 EVERY TIME. */
322 : /* MISSLX(J) = MALLOW EXCEPT WHEN A GROUP IS ALL ONE VALUE (AND */
323 : /* LBIT(J) = 0) AND THAT VALUE IS MISSING. IN */
324 : /* THAT CASE, MISSLX(J) IS MISSP OR MISSS. THIS */
325 : /* GETS INSERTED INTO JMIN(J) LATER AS THE */
326 : /* MISSING INDICATOR; IT CAN'T BE PUT IN UNTIL */
327 : /* THE END, BECAUSE JMIN( ) IS USED TO CALCULATE */
328 : /* THE MAXIMUM NUMBER OF BITS (IBITS) NEEDED TO */
329 : /* PACK JMIN( ). */
330 : /* 1 2 3 4 5 6 7 X */
331 :
332 : /* NON SYSTEM SUBROUTINES CALLED */
333 : /* NONE */
334 :
335 :
336 :
337 : /* MISSLX( ) was AN AUTOMATIC ARRAY. */
338 13 : misslx = (integer *)calloc(*ndg,sizeof(integer));
339 13 : if( misslx == NULL )
340 : {
341 0 : *ier = -1;
342 0 : return 0;
343 : }
344 :
345 : /* Parameter adjustments */
346 13 : --ic;
347 13 : --nov;
348 13 : --lbit;
349 13 : --jmax;
350 13 : --jmin;
351 :
352 : /* Function Body */
353 :
354 13 : *ier = 0;
355 13 : iersav = 0;
356 : /* CALL TIMPR(KFILDO,KFILDO,'START PACK_GP ') */
357 13 : *(unsigned char *)cfeed = (char) ifeed;
358 :
359 13 : ired = 0;
360 : /* IRED IS A FLAG. WHEN ZERO, REDUCE WILL BE CALLED. */
361 : /* IF REDUCE ABORTS, IRED = 1 AND IS NOT CALLED. IN */
362 : /* THIS CASE PACK_GP EXECUTES AGAIN EXCEPT FOR REDUCE. */
363 :
364 13 : if (*inc <= 0) {
365 0 : iersav = 717;
366 : /* WRITE(KFILDO,101)INC */
367 : /* 101 FORMAT(/' ****INC ='I8,' NOT CORRECT IN PACK_GP. 1 IS USED.') */
368 : }
369 :
370 : /* THERE WILL BE A RESTART OF PACK_GP IF SUBROUTINE REDUCE */
371 : /* ABORTS. THIS SHOULD NOT HAPPEN, BUT IF IT DOES, PACK_GP */
372 : /* WILL COMPLETE WITHOUT SUBROUTINE REDUCE. A NON FATAL */
373 : /* DIAGNOSTIC RETURN IS PROVIDED. */
374 :
375 13 : L102:
376 : /*kinc = max(*inc,1);*/
377 14 : kinc = (*inc > 1) ? *inc : 1;
378 14 : lminpk = *minpk;
379 :
380 : /* CALCULATE THE POWERS OF 2 THE FIRST TIME ENTERED. */
381 :
382 14 : if (ifirst == 0) {
383 13 : ifirst = 1;
384 13 : ibxx2[0] = 1;
385 :
386 403 : for (j = 1; j <= 30; ++j) {
387 390 : ibxx2[j] = ibxx2[j - 1] << 1;
388 : /* L104: */
389 : }
390 :
391 : }
392 :
393 : /* THERE WILL BE A RESTART AT 105 IS NDG IS NOT LARGE ENOUGH. */
394 : /* A NON FATAL DIAGNOSTIC RETURN IS PROVIDED. */
395 :
396 14 : L105:
397 14 : kstart = 1;
398 14 : ktotal = 0;
399 14 : *lx = 0;
400 14 : adda = FALSE_;
401 14 : lmiss = 0;
402 14 : if (*is523 == 1) {
403 4 : lmiss = 1;
404 : }
405 14 : if (*is523 == 2) {
406 0 : lmiss = 2;
407 : }
408 :
409 : /* ************************************* */
410 :
411 : /* THIS SECTION COMPUTES STATISTICS FOR GROUP A. GROUP A IS */
412 : /* A GROUP OF SIZE LMINPK. */
413 :
414 : /* ************************************* */
415 :
416 14 : ibita = 0;
417 14 : mina = mallow;
418 14 : maxa = -mallow;
419 14 : minak = mallow;
420 14 : maxak = -mallow;
421 :
422 : /* FIND THE MIN AND MAX OF GROUP A. THIS WILL INITIALLY BE OF */
423 : /* SIZE LMINPK (IF THERE ARE STILL LMINPK VALUES IN IC( )), BUT */
424 : /* WILL INCREASE IN SIZE IN INCREMENTS OF INC UNTIL A NEW */
425 : /* GROUP IS STARTED. THE DEFINITION OF GROUP A IS DONE HERE */
426 : /* ONLY ONCE (UPON INITIAL ENTRY), BECAUSE A GROUP B CAN ALWAYS */
427 : /* BECOME A NEW GROUP A AFTER A IS PACKED, EXCEPT IF LMINPK */
428 : /* HAS TO BE INCREASED BECAUSE NDG IS TOO SMALL. THEREFORE, */
429 : /* THE SEPARATE LOOPS FOR MISSING AND NON-MISSING HERE BUYS */
430 : /* ALMOST NOTHING. */
431 :
432 : /* Computing MIN */
433 14 : i__1 = kstart + lminpk - 1;
434 : /*nenda = min(i__1,*nxy);*/
435 14 : nenda = (i__1 < *nxy) ? i__1 : *nxy;
436 14 : if (*nxy - nenda <= lminpk / 2) {
437 2 : nenda = *nxy;
438 : }
439 : /* ABOVE STATEMENT GUARANTEES THE LAST GROUP IS GT LMINPK/2 BY */
440 : /* MAKING THE ACTUAL GROUP LARGER. IF A PROVISION LIKE THIS IS */
441 : /* NOT INCLUDED, THERE WILL MANY TIMES BE A VERY SMALL GROUP */
442 : /* AT THE END. USE SEPARATE LOOPS FOR MISSING AND NO MISSING */
443 : /* VALUES FOR EFFICIENCY. */
444 :
445 : /* DETERMINE WHETHER THERE IS A LONG STRING OF THE SAME VALUE */
446 : /* UNLESS NENDA = NXY. THIS MAY ALLOW A LARGE GROUP A TO */
447 : /* START WITH, AS WITH MISSING VALUES. SEPARATE LOOPS FOR */
448 : /* MISSING OPTIONS. THIS SECTION IS ONLY EXECUTED ONCE, */
449 : /* IN DETERMINING THE FIRST GROUP. IT HELPS FOR AN ARRAY */
450 : /* OF MOSTLY MISSING VALUES OR OF ONE VALUE, SUCH AS */
451 : /* RADAR OR PRECIP DATA. */
452 :
453 14 : if (nenda != *nxy && ic[kstart] == ic[kstart + 1]) {
454 : /* NO NEED TO EXECUTE IF FIRST TWO VALUES ARE NOT EQUAL. */
455 :
456 10 : if (*is523 == 0) {
457 : /* THIS LOOP IS FOR NO MISSING VALUES. */
458 :
459 6 : i__1 = *nxy;
460 12 : for (k = kstart + 1; k <= i__1; ++k) {
461 :
462 12 : if (ic[k] != ic[kstart]) {
463 : /* Computing MAX */
464 6 : i__2 = nenda;
465 6 : i__3 = k - 1;
466 : /*nenda = max(i__2,i__3);*/
467 6 : nenda = (i__2 > i__3) ? i__2 : i__3;
468 6 : goto L114;
469 : }
470 :
471 : /* L111: */
472 : }
473 :
474 0 : nenda = *nxy;
475 : /* FALL THROUGH THE LOOP MEANS ALL VALUES ARE THE SAME. */
476 :
477 4 : } else if (*is523 == 1) {
478 : /* THIS LOOP IS FOR PRIMARY MISSING VALUES ONLY. */
479 :
480 4 : i__1 = *nxy;
481 82580 : for (k = kstart + 1; k <= i__1; ++k) {
482 :
483 82580 : if (ic[k] != *missp) {
484 :
485 4 : if (ic[k] != ic[kstart]) {
486 : /* Computing MAX */
487 4 : i__2 = nenda;
488 4 : i__3 = k - 1;
489 : /*nenda = max(i__2,i__3);*/
490 4 : nenda = (i__2 > i__3) ? i__2 : i__3;
491 4 : goto L114;
492 : }
493 :
494 : }
495 :
496 : /* L112: */
497 : }
498 :
499 0 : nenda = *nxy;
500 : /* FALL THROUGH THE LOOP MEANS ALL VALUES ARE THE SAME. */
501 :
502 : } else {
503 : /* THIS LOOP IS FOR PRIMARY AND SECONDARY MISSING VALUES. */
504 :
505 0 : i__1 = *nxy;
506 0 : for (k = kstart + 1; k <= i__1; ++k) {
507 :
508 0 : if (ic[k] != *missp && ic[k] != *misss) {
509 :
510 0 : if (ic[k] != ic[kstart]) {
511 : /* Computing MAX */
512 0 : i__2 = nenda;
513 0 : i__3 = k - 1;
514 : /*nenda = max(i__2,i__3);*/
515 0 : nenda = (i__2 > i__3) ? i__2 : i__3;
516 0 : goto L114;
517 : }
518 :
519 : }
520 :
521 : /* L113: */
522 : }
523 :
524 0 : nenda = *nxy;
525 : /* FALL THROUGH THE LOOP MEANS ALL VALUES ARE THE SAME. */
526 : }
527 :
528 : }
529 :
530 4 : L114:
531 14 : if (*is523 == 0) {
532 :
533 10 : i__1 = nenda;
534 98 : for (k = kstart; k <= i__1; ++k) {
535 88 : if (ic[k] < mina) {
536 27 : mina = ic[k];
537 27 : minak = k;
538 : }
539 88 : if (ic[k] > maxa) {
540 22 : maxa = ic[k];
541 22 : maxak = k;
542 : }
543 : /* L115: */
544 : }
545 :
546 4 : } else if (*is523 == 1) {
547 :
548 4 : i__1 = nenda;
549 82584 : for (k = kstart; k <= i__1; ++k) {
550 82580 : if (ic[k] == *missp) {
551 82580 : goto L117;
552 : }
553 0 : if (ic[k] < mina) {
554 0 : mina = ic[k];
555 0 : minak = k;
556 : }
557 0 : if (ic[k] > maxa) {
558 0 : maxa = ic[k];
559 0 : maxak = k;
560 : }
561 82580 : L117:
562 : ;
563 : }
564 :
565 : } else {
566 :
567 0 : i__1 = nenda;
568 0 : for (k = kstart; k <= i__1; ++k) {
569 0 : if (ic[k] == *missp || ic[k] == *misss) {
570 0 : goto L120;
571 : }
572 0 : if (ic[k] < mina) {
573 0 : mina = ic[k];
574 0 : minak = k;
575 : }
576 0 : if (ic[k] > maxa) {
577 0 : maxa = ic[k];
578 0 : maxak = k;
579 : }
580 0 : L120:
581 : ;
582 : }
583 :
584 : }
585 :
586 14 : kounta = nenda - kstart + 1;
587 :
588 : /* INCREMENT KTOTAL AND FIND THE BITS NEEDED TO PACK THE A GROUP. */
589 :
590 14 : ktotal += kounta;
591 14 : mislla = 0;
592 14 : if (mina != mallow) {
593 10 : goto L125;
594 : }
595 : /* ALL MISSING VALUES MUST BE ACCOMMODATED. */
596 4 : mina = 0;
597 4 : maxa = 0;
598 4 : mislla = 1;
599 4 : ibitb = 0;
600 4 : if (*is523 != 2) {
601 4 : goto L130;
602 : }
603 : /* WHEN ALL VALUES ARE MISSING AND THERE ARE NO */
604 : /* SECONDARY MISSING VALUES, IBITA = 0. */
605 : /* OTHERWISE, IBITA MUST BE CALCULATED. */
606 :
607 0 : L125:
608 10 : itest = maxa - mina + lmiss;
609 :
610 74 : for (ibita = 0; ibita <= 30; ++ibita) {
611 74 : if (itest < ibxx2[ibita]) {
612 10 : goto L130;
613 : }
614 : /* *** THIS TEST IS THE SAME AS: */
615 : /* *** IF(MAXA-MINA.LT.IBXX2(IBITA)-LMISS)GO TO 130 */
616 : /* L126: */
617 : }
618 :
619 : /* WRITE(KFILDO,127)MAXA,MINA */
620 : /* 127 FORMAT(' ****ERROR IN PACK_GP. VALUE WILL NOT PACK IN 30 BITS.', */
621 : /* 1 ' MAXA ='I13,' MINA ='I13,'. ERROR AT 127.') */
622 0 : *ier = 706;
623 0 : goto L900;
624 :
625 39979 : L130:
626 :
627 : /* ***D WRITE(KFILDO,131)KOUNTA,KTOTAL,MINA,MAXA,IBITA,MISLLA */
628 : /* ***D131 FORMAT(' AT 130, KOUNTA ='I8,' KTOTAL ='I8,' MINA ='I8, */
629 : /* ***D 1 ' MAXA ='I8,' IBITA ='I3,' MISLLA ='I3) */
630 :
631 39979 : L133:
632 39979 : if (ktotal >= *nxy) {
633 5 : goto L200;
634 : }
635 :
636 : /* ************************************* */
637 :
638 : /* THIS SECTION COMPUTES STATISTICS FOR GROUP B. GROUP B IS A */
639 : /* GROUP OF SIZE LMINPK IMMEDIATELY FOLLOWING GROUP A. */
640 :
641 : /* ************************************* */
642 :
643 39974 : L140:
644 526217 : minb = mallow;
645 526217 : maxb = -mallow;
646 526217 : minbk = mallow;
647 526217 : maxbk = -mallow;
648 526217 : ibitbs = 0;
649 526217 : mstart = ktotal + 1;
650 :
651 : /* DETERMINE WHETHER THERE IS A LONG STRING OF THE SAME VALUE. */
652 : /* THIS WORKS WHEN THERE ARE NO MISSING VALUES. */
653 :
654 526217 : nendb = 1;
655 :
656 526217 : if (mstart < *nxy) {
657 :
658 526217 : if (*is523 == 0) {
659 : /* THIS LOOP IS FOR NO MISSING VALUES. */
660 :
661 344 : i__1 = *nxy;
662 379 : for (k = mstart + 1; k <= i__1; ++k) {
663 :
664 379 : if (ic[k] != ic[mstart]) {
665 344 : nendb = k - 1;
666 344 : goto L150;
667 : }
668 :
669 : /* L145: */
670 : }
671 :
672 0 : nendb = *nxy;
673 : /* FALL THROUGH THE LOOP MEANS ALL REMAINING VALUES */
674 : /* ARE THE SAME. */
675 : }
676 :
677 : }
678 :
679 525873 : L150:
680 : /* Computing MAX */
681 : /* Computing MIN */
682 645128 : i__3 = ktotal + lminpk;
683 : /*i__1 = nendb, i__2 = min(i__3,*nxy);*/
684 645128 : i__1 = nendb;
685 645128 : i__2 = (i__3 < *nxy) ? i__3 : *nxy;
686 : /*nendb = max(i__1,i__2);*/
687 645128 : nendb = (i__1 > i__2) ? i__1 : i__2;
688 : /* **** 150 NENDB=MIN(KTOTAL+LMINPK,NXY) */
689 :
690 645128 : if (*nxy - nendb <= lminpk / 2) {
691 93 : nendb = *nxy;
692 : }
693 : /* ABOVE STATEMENT GUARANTEES THE LAST GROUP IS GT LMINPK/2 BY */
694 : /* MAKING THE ACTUAL GROUP LARGER. IF A PROVISION LIKE THIS IS */
695 : /* NOT INCLUDED, THERE WILL MANY TIMES BE A VERY SMALL GROUP */
696 : /* AT THE END. USE SEPARATE LOOPS FOR MISSING AND NO MISSING */
697 :
698 : /* USE SEPARATE LOOPS FOR MISSING AND NO MISSING VALUES */
699 : /* FOR EFFICIENCY. */
700 :
701 645128 : if (*is523 == 0) {
702 :
703 1262 : i__1 = nendb;
704 5583 : for (k = mstart; k <= i__1; ++k) {
705 4321 : if (ic[k] <= minb) {
706 1336 : minb = ic[k];
707 : /* NOTE LE, NOT LT. LT COULD BE USED BUT THEN A */
708 : /* RECOMPUTE OVER THE WHOLE GROUP WOULD BE NEEDED */
709 : /* MORE OFTEN. SAME REASONING FOR GE AND OTHER */
710 : /* LOOPS BELOW. */
711 1336 : minbk = k;
712 : }
713 4321 : if (ic[k] >= maxb) {
714 1146 : maxb = ic[k];
715 1146 : maxbk = k;
716 : }
717 : /* L155: */
718 : }
719 :
720 643866 : } else if (*is523 == 1) {
721 :
722 643866 : i__1 = nendb;
723 6020600 : for (k = mstart; k <= i__1; ++k) {
724 5376730 : if (ic[k] == *missp) {
725 4740240 : goto L157;
726 : }
727 636492 : if (ic[k] <= minb) {
728 300089 : minb = ic[k];
729 300089 : minbk = k;
730 : }
731 636492 : if (ic[k] >= maxb) {
732 309445 : maxb = ic[k];
733 309445 : maxbk = k;
734 : }
735 5376730 : L157:
736 : ;
737 : }
738 :
739 : } else {
740 :
741 0 : i__1 = nendb;
742 0 : for (k = mstart; k <= i__1; ++k) {
743 0 : if (ic[k] == *missp || ic[k] == *misss) {
744 0 : goto L160;
745 : }
746 0 : if (ic[k] <= minb) {
747 0 : minb = ic[k];
748 0 : minbk = k;
749 : }
750 0 : if (ic[k] >= maxb) {
751 0 : maxb = ic[k];
752 0 : maxbk = k;
753 : }
754 0 : L160:
755 : ;
756 : }
757 :
758 : }
759 :
760 645128 : kountb = nendb - ktotal;
761 645128 : misllb = 0;
762 645128 : if (minb != mallow) {
763 181943 : goto L165;
764 : }
765 : /* ALL MISSING VALUES MUST BE ACCOMMODATED. */
766 463185 : minb = 0;
767 463185 : maxb = 0;
768 463185 : misllb = 1;
769 463185 : ibitb = 0;
770 :
771 463185 : if (*is523 != 2) {
772 463185 : goto L170;
773 : }
774 : /* WHEN ALL VALUES ARE MISSING AND THERE ARE NO SECONDARY */
775 : /* MISSING VALUES, IBITB = 0. OTHERWISE, IBITB MUST BE */
776 : /* CALCULATED. */
777 :
778 0 : L165:
779 181943 : if( (GIntBig)maxb - minb < INT_MIN ||
780 181943 : (GIntBig)maxb - minb > INT_MAX )
781 : {
782 0 : *ier = -1;
783 0 : free(misslx);
784 0 : return 0;
785 : }
786 :
787 439657 : for (ibitb = ibitbs; ibitb <= 30; ++ibitb) {
788 439657 : if (maxb - minb < ibxx2[ibitb] - lmiss) {
789 181943 : goto L170;
790 : }
791 : /* L166: */
792 : }
793 :
794 : /* WRITE(KFILDO,167)MAXB,MINB */
795 : /* 167 FORMAT(' ****ERROR IN PACK_GP. VALUE WILL NOT PACK IN 30 BITS.', */
796 : /* 1 ' MAXB ='I13,' MINB ='I13,'. ERROR AT 167.') */
797 0 : *ier = 706;
798 0 : goto L900;
799 :
800 : /* COMPARE THE BITS NEEDED TO PACK GROUP B WITH THOSE NEEDED */
801 : /* TO PACK GROUP A. IF IBITB GE IBITA, TRY TO ADD TO GROUP A. */
802 : /* IF NOT, TRY TO ADD A'S POINTS TO B, UNLESS ADDITION TO A */
803 : /* HAS BEEN DONE. THIS LATTER IS CONTROLLED WITH ADDA. */
804 :
805 645128 : L170:
806 :
807 : /* ***D WRITE(KFILDO,171)KOUNTA,KTOTAL,MINA,MAXA,IBITA,MISLLA, */
808 : /* ***D 1 MINB,MAXB,IBITB,MISLLB */
809 : /* ***D171 FORMAT(' AT 171, KOUNTA ='I8,' KTOTAL ='I8,' MINA ='I8, */
810 : /* ***D 1 ' MAXA ='I8,' IBITA ='I3,' MISLLA ='I3, */
811 : /* ***D 2 ' MINB ='I8,' MAXB ='I8,' IBITB ='I3,' MISLLB ='I3) */
812 :
813 645128 : if (ibitb >= ibita) {
814 625900 : goto L180;
815 : }
816 19228 : if (adda) {
817 4867 : goto L200;
818 : }
819 :
820 : /* ************************************* */
821 :
822 : /* GROUP B REQUIRES LESS BITS THAN GROUP A. PUT AS MANY OF A'S */
823 : /* POINTS INTO B AS POSSIBLE WITHOUT EXCEEDING THE NUMBER OF */
824 : /* BITS NECESSARY TO PACK GROUP B. */
825 :
826 : /* ************************************* */
827 :
828 14361 : kounts = kounta;
829 : /* KOUNTA REFERS TO THE PRESENT GROUP A. */
830 14361 : mintst = minb;
831 14361 : maxtst = maxb;
832 14361 : mintstk = minbk;
833 14361 : maxtstk = maxbk;
834 :
835 : /* USE SEPARATE LOOPS FOR MISSING AND NO MISSING VALUES */
836 : /* FOR EFFICIENCY. */
837 :
838 14361 : if (*is523 == 0) {
839 :
840 57 : i__1 = kstart;
841 417 : for (k = ktotal; k >= i__1; --k) {
842 : /* START WITH THE END OF THE GROUP AND WORK BACKWARDS. */
843 417 : if (ic[k] < minb) {
844 24 : mintst = ic[k];
845 24 : mintstk = k;
846 393 : } else if (ic[k] > maxb) {
847 55 : maxtst = ic[k];
848 55 : maxtstk = k;
849 : }
850 417 : if (maxtst - mintst >= ibxx2[ibitb]) {
851 57 : goto L174;
852 : }
853 : /* NOTE THAT FOR THIS LOOP, LMISS = 0. */
854 360 : minb = mintst;
855 360 : maxb = maxtst;
856 360 : minbk = mintstk;
857 360 : maxbk = maxtstk;
858 360 : --kounta;
859 : /* THERE IS ONE LESS POINT NOW IN A. */
860 : /* L1715: */
861 : }
862 :
863 14304 : } else if (*is523 == 1) {
864 :
865 14304 : i__1 = kstart;
866 68552 : for (k = ktotal; k >= i__1; --k) {
867 : /* START WITH THE END OF THE GROUP AND WORK BACKWARDS. */
868 68552 : if (ic[k] == *missp) {
869 25363 : goto L1718;
870 : }
871 43189 : if (ic[k] < minb) {
872 7648 : mintst = ic[k];
873 7648 : mintstk = k;
874 35541 : } else if (ic[k] > maxb) {
875 12261 : maxtst = ic[k];
876 12261 : maxtstk = k;
877 : }
878 43189 : if (maxtst - mintst >= ibxx2[ibitb] - lmiss) {
879 14304 : goto L174;
880 : }
881 : /* FOR THIS LOOP, LMISS = 1. */
882 28885 : minb = mintst;
883 28885 : maxb = maxtst;
884 28885 : minbk = mintstk;
885 28885 : maxbk = maxtstk;
886 28885 : misllb = 0;
887 : /* WHEN THE POINT IS NON MISSING, MISLLB SET = 0. */
888 54248 : L1718:
889 54248 : --kounta;
890 : /* THERE IS ONE LESS POINT NOW IN A. */
891 : /* L1719: */
892 : }
893 :
894 : } else {
895 :
896 0 : i__1 = kstart;
897 0 : for (k = ktotal; k >= i__1; --k) {
898 : /* START WITH THE END OF THE GROUP AND WORK BACKWARDS. */
899 0 : if (ic[k] == *missp || ic[k] == *misss) {
900 0 : goto L1729;
901 : }
902 0 : if (ic[k] < minb) {
903 0 : mintst = ic[k];
904 0 : mintstk = k;
905 0 : } else if (ic[k] > maxb) {
906 0 : maxtst = ic[k];
907 0 : maxtstk = k;
908 : }
909 0 : if (maxtst - mintst >= ibxx2[ibitb] - lmiss) {
910 0 : goto L174;
911 : }
912 : /* FOR THIS LOOP, LMISS = 2. */
913 0 : minb = mintst;
914 0 : maxb = maxtst;
915 0 : minbk = mintstk;
916 0 : maxbk = maxtstk;
917 0 : misllb = 0;
918 : /* WHEN THE POINT IS NON MISSING, MISLLB SET = 0. */
919 0 : L1729:
920 0 : --kounta;
921 : /* THERE IS ONE LESS POINT NOW IN A. */
922 : /* L173: */
923 : }
924 :
925 : }
926 :
927 : /* AT THIS POINT, KOUNTA CONTAINS THE NUMBER OF POINTS TO CLOSE */
928 : /* OUT GROUP A WITH. GROUP B NOW STARTS WITH KSTART+KOUNTA AND */
929 : /* ENDS WITH NENDB. MINB AND MAXB HAVE BEEN ADJUSTED AS */
930 : /* NECESSARY TO REFLECT GROUP B (EVEN THOUGH THE NUMBER OF BITS */
931 : /* NEEDED TO PACK GROUP B HAVE NOT INCREASED, THE END POINTS */
932 : /* OF THE RANGE MAY HAVE). */
933 :
934 0 : L174:
935 14361 : if (kounta == kounts) {
936 2286 : goto L200;
937 : }
938 : /* ON TRANSFER, GROUP A WAS NOT CHANGED. CLOSE IT OUT. */
939 :
940 : /* ONE OR MORE POINTS WERE TAKEN OUT OF A. RANGE AND IBITA */
941 : /* MAY HAVE TO BE RECOMPUTED; IBITA COULD BE LESS THAN */
942 : /* ORIGINALLY COMPUTED. IN FACT, GROUP A CAN NOW CONTAIN */
943 : /* ONLY ONE POINT AND BE PACKED WITH ZERO BITS */
944 : /* (UNLESS MISSS NE 0). */
945 :
946 12075 : nouta = kounts - kounta;
947 12075 : ktotal -= nouta;
948 12075 : kountb += nouta;
949 12075 : if (nenda - nouta > minak && nenda - nouta > maxak) {
950 1850 : goto L200;
951 : }
952 : /* WHEN THE ABOVE TEST IS MET, THE MIN AND MAX OF THE */
953 : /* CURRENT GROUP A WERE WITHIN THE OLD GROUP A, SO THE */
954 : /* RANGE AND IBITA DO NOT NEED TO BE RECOMPUTED. */
955 : /* NOTE THAT MINAK AND MAXAK ARE NO LONGER NEEDED. */
956 10225 : ibita = 0;
957 10225 : mina = mallow;
958 10225 : maxa = -mallow;
959 :
960 : /* USE SEPARATE LOOPS FOR MISSING AND NO MISSING VALUES */
961 : /* FOR EFFICIENCY. */
962 :
963 10225 : if (*is523 == 0) {
964 :
965 51 : i__1 = nenda - nouta;
966 246 : for (k = kstart; k <= i__1; ++k) {
967 195 : if (ic[k] < mina) {
968 92 : mina = ic[k];
969 : }
970 195 : if (ic[k] > maxa) {
971 89 : maxa = ic[k];
972 : }
973 : /* L1742: */
974 : }
975 :
976 10174 : } else if (*is523 == 1) {
977 :
978 10174 : i__1 = nenda - nouta;
979 69551 : for (k = kstart; k <= i__1; ++k) {
980 59377 : if (ic[k] == *missp) {
981 5460 : goto L1744;
982 : }
983 53917 : if (ic[k] < mina) {
984 21942 : mina = ic[k];
985 : }
986 53917 : if (ic[k] > maxa) {
987 21347 : maxa = ic[k];
988 : }
989 59377 : L1744:
990 : ;
991 : }
992 :
993 : } else {
994 :
995 0 : i__1 = nenda - nouta;
996 0 : for (k = kstart; k <= i__1; ++k) {
997 0 : if (ic[k] == *missp || ic[k] == *misss) {
998 0 : goto L175;
999 : }
1000 0 : if (ic[k] < mina) {
1001 0 : mina = ic[k];
1002 : }
1003 0 : if (ic[k] > maxa) {
1004 0 : maxa = ic[k];
1005 : }
1006 0 : L175:
1007 : ;
1008 : }
1009 :
1010 : }
1011 :
1012 10225 : mislla = 0;
1013 10225 : if (mina != mallow) {
1014 10225 : goto L1750;
1015 : }
1016 : /* ALL MISSING VALUES MUST BE ACCOMMODATED. */
1017 0 : mina = 0;
1018 0 : maxa = 0;
1019 0 : mislla = 1;
1020 0 : if (*is523 != 2) {
1021 0 : goto L177;
1022 : }
1023 : /* WHEN ALL VALUES ARE MISSING AND THERE ARE NO SECONDARY */
1024 : /* MISSING VALUES IBITA = 0 AS ORIGINALLY SET. OTHERWISE, */
1025 : /* IBITA MUST BE CALCULATED. */
1026 :
1027 0 : L1750:
1028 10225 : itest = maxa - mina + lmiss;
1029 :
1030 44270 : for (ibita = 0; ibita <= 30; ++ibita) {
1031 44270 : if (itest < ibxx2[ibita]) {
1032 10225 : goto L177;
1033 : }
1034 : /* *** THIS TEST IS THE SAME AS: */
1035 : /* *** IF(MAXA-MINA.LT.IBXX2(IBITA)-LMISS)GO TO 177 */
1036 : /* L176: */
1037 : }
1038 :
1039 : /* WRITE(KFILDO,1760)MAXA,MINA */
1040 : /* 1760 FORMAT(' ****ERROR IN PACK_GP. VALUE WILL NOT PACK IN 30 BITS.', */
1041 : /* 1 ' MAXA ='I13,' MINA ='I13,'. ERROR AT 1760.') */
1042 0 : *ier = 706;
1043 0 : goto L900;
1044 :
1045 10225 : L177:
1046 10225 : goto L200;
1047 :
1048 : /* ************************************* */
1049 :
1050 : /* AT THIS POINT, GROUP B REQUIRES AS MANY BITS TO PACK AS GROUPA. */
1051 : /* THEREFORE, TRY TO ADD INC POINTS TO GROUP A WITHOUT INCREASING */
1052 : /* IBITA. THIS AUGMENTED GROUP IS CALLED GROUP C. */
1053 :
1054 : /* ************************************* */
1055 :
1056 625900 : L180:
1057 625900 : if (mislla == 1) {
1058 492551 : minc = mallow;
1059 492551 : minck = mallow;
1060 492551 : maxc = -mallow;
1061 492551 : maxck = -mallow;
1062 : } else {
1063 133349 : minc = mina;
1064 133349 : maxc = maxa;
1065 133349 : minck = minak;
1066 133349 : maxck = minak;
1067 : }
1068 :
1069 625900 : nount = 0;
1070 625900 : if (*nxy - (ktotal + kinc) <= lminpk / 2) {
1071 12 : kinc = *nxy - ktotal;
1072 : }
1073 : /* ABOVE STATEMENT CONSTRAINS THE LAST GROUP TO BE NOT LESS THAN */
1074 : /* LMINPK/2 IN SIZE. IF A PROVISION LIKE THIS IS NOT INCLUDED, */
1075 : /* THERE WILL MANY TIMES BE A VERY SMALL GROUP AT THE END. */
1076 :
1077 : /* USE SEPARATE LOOPS FOR MISSING AND NO MISSING VALUES */
1078 : /* FOR EFFICIENCY. SINCE KINC IS USUALLY 1, USING SEPARATE */
1079 : /* LOOPS HERE DOESN'T BUY MUCH. A MISSING VALUE WILL ALWAYS */
1080 : /* TRANSFER BACK TO GROUP A. */
1081 :
1082 625900 : if (*is523 == 0) {
1083 :
1084 : /* Computing MIN */
1085 1163 : i__2 = ktotal + kinc;
1086 : /*i__1 = min(i__2,*nxy);*/
1087 1163 : i__1 = (i__2 < *nxy) ? i__2 : *nxy;
1088 2366 : for (k = ktotal + 1; k <= i__1; ++k) {
1089 1203 : if (ic[k] < minc) {
1090 59 : minc = ic[k];
1091 59 : minck = k;
1092 : }
1093 1203 : if (ic[k] > maxc) {
1094 99 : maxc = ic[k];
1095 99 : maxck = k;
1096 : }
1097 1203 : ++nount;
1098 : /* L185: */
1099 : }
1100 :
1101 624737 : } else if (*is523 == 1) {
1102 :
1103 : /* Computing MIN */
1104 624737 : i__2 = ktotal + kinc;
1105 : /*i__1 = min(i__2,*nxy);*/
1106 624737 : i__1 = (i__2 < *nxy) ? i__2 : *nxy;
1107 1249490 : for (k = ktotal + 1; k <= i__1; ++k) {
1108 624757 : if (ic[k] == *missp) {
1109 498577 : goto L186;
1110 : }
1111 126180 : if (ic[k] < minc) {
1112 16225 : minc = ic[k];
1113 16225 : minck = k;
1114 : }
1115 126180 : if (ic[k] > maxc) {
1116 19974 : maxc = ic[k];
1117 19974 : maxck = k;
1118 : }
1119 106206 : L186:
1120 624757 : ++nount;
1121 : /* L187: */
1122 : }
1123 :
1124 : } else {
1125 :
1126 : /* Computing MIN */
1127 0 : i__2 = ktotal + kinc;
1128 : /*i__1 = min(i__2,*nxy);*/
1129 0 : i__1 = (i__2 < *nxy) ? i__2 : *nxy;
1130 0 : for (k = ktotal + 1; k <= i__1; ++k) {
1131 0 : if (ic[k] == *missp || ic[k] == *misss) {
1132 0 : goto L189;
1133 : }
1134 0 : if (ic[k] < minc) {
1135 0 : minc = ic[k];
1136 0 : minck = k;
1137 : }
1138 0 : if (ic[k] > maxc) {
1139 0 : maxc = ic[k];
1140 0 : maxck = k;
1141 : }
1142 0 : L189:
1143 0 : ++nount;
1144 : /* L190: */
1145 : }
1146 :
1147 : }
1148 :
1149 : /* ***D WRITE(KFILDO,191)KOUNTA,KTOTAL,MINA,MAXA,IBITA,MISLLA, */
1150 : /* ***D 1 MINC,MAXC,NOUNT,IC(KTOTAL),IC(KTOTAL+1) */
1151 : /* ***D191 FORMAT(' AT 191, KOUNTA ='I8,' KTOTAL ='I8,' MINA ='I8, */
1152 : /* ***D 1 ' MAXA ='I8,' IBITA ='I3,' MISLLA ='I3, */
1153 : /* ***D 2 ' MINC ='I8,' MAXC ='I8, */
1154 : /* ***D 3 ' NOUNT ='I5,' IC(KTOTAL) ='I9,' IC(KTOTAL+1) =',I9) */
1155 :
1156 : /* IF THE NUMBER OF BITS NEEDED FOR GROUP C IS GT IBITA, */
1157 : /* THEN THIS GROUP A IS A GROUP TO PACK. */
1158 :
1159 625900 : if (minc == mallow) {
1160 488680 : minc = mina;
1161 488680 : maxc = maxa;
1162 488680 : minck = minak;
1163 488680 : maxck = maxak;
1164 488680 : misllc = 1;
1165 488680 : goto L195;
1166 : /* WHEN THE NEW VALUE(S) ARE MISSING, THEY CAN ALWAYS */
1167 : /* BE ADDED. */
1168 :
1169 : } else {
1170 137220 : misllc = 0;
1171 : }
1172 :
1173 137220 : if (maxc - minc >= ibxx2[ibita] - lmiss) {
1174 20737 : goto L200;
1175 : }
1176 :
1177 : /* THE BITS NECESSARY FOR GROUP C HAS NOT INCREASED FROM THE */
1178 : /* BITS NECESSARY FOR GROUP A. ADD THIS POINT(S) TO GROUP A. */
1179 : /* COMPUTE THE NEXT GROUP B, ETC., UNLESS ALL POINTS HAVE BEEN */
1180 : /* USED. */
1181 :
1182 116483 : L195:
1183 605163 : ktotal += nount;
1184 605163 : kounta += nount;
1185 605163 : mina = minc;
1186 605163 : maxa = maxc;
1187 605163 : minak = minck;
1188 605163 : maxak = maxck;
1189 605163 : mislla = misllc;
1190 605163 : adda = TRUE_;
1191 605163 : if (ktotal >= *nxy) {
1192 9 : goto L200;
1193 : }
1194 :
1195 605154 : if (minbk > ktotal && maxbk > ktotal) {
1196 118911 : mstart = nendb + 1;
1197 : /* THE MAX AND MIN OF GROUP B WERE NOT FROM THE POINTS */
1198 : /* REMOVED, SO THE WHOLE GROUP DOES NOT HAVE TO BE LOOKED */
1199 : /* AT TO DETERMINE THE NEW MAX AND MIN. RATHER START */
1200 : /* JUST BEYOND THE OLD NENDB. */
1201 118911 : ibitbs = ibitb;
1202 118911 : nendb = 1;
1203 118911 : goto L150;
1204 : } else {
1205 486243 : goto L140;
1206 : }
1207 :
1208 : /* ************************************* */
1209 :
1210 : /* GROUP A IS TO BE PACKED. STORE VALUES IN JMIN( ), JMAX( ), */
1211 : /* LBIT( ), AND NOV( ). */
1212 :
1213 : /* ************************************* */
1214 :
1215 39979 : L200:
1216 39979 : ++(*lx);
1217 39979 : if (*lx <= *ndg) {
1218 39979 : goto L205;
1219 : }
1220 0 : lminpk += lminpk / 2;
1221 : /* WRITE(KFILDO,201)NDG,LMINPK,LX */
1222 : /* 201 FORMAT(' ****NDG ='I5,' NOT LARGE ENOUGH.', */
1223 : /* 1 ' LMINPK IS INCREASED TO 'I3,' FOR THIS FIELD.'/ */
1224 : /* 2 ' LX = 'I10) */
1225 0 : iersav = 716;
1226 0 : goto L105;
1227 :
1228 39979 : L205:
1229 39979 : jmin[*lx] = mina;
1230 39979 : jmax[*lx] = maxa;
1231 39979 : lbit[*lx] = ibita;
1232 39979 : nov[*lx] = kounta;
1233 39979 : kstart = ktotal + 1;
1234 :
1235 39979 : if (mislla == 0) {
1236 36105 : misslx[*lx - 1] = mallow;
1237 : } else {
1238 3874 : misslx[*lx - 1] = ic[ktotal];
1239 : /* IC(KTOTAL) WAS THE LAST VALUE PROCESSED. IF MISLLA NE 0, */
1240 : /* THIS MUST BE THE MISSING VALUE FOR THIS GROUP. */
1241 : }
1242 :
1243 : /* ***D WRITE(KFILDO,206)MISLLA,IC(KTOTAL),KTOTAL,LX,JMIN(LX),JMAX(LX), */
1244 : /* ***D 1 LBIT(LX),NOV(LX),MISSLX(LX) */
1245 : /* ***D206 FORMAT(' AT 206, MISLLA ='I2,' IC(KTOTAL) ='I5,' KTOTAL ='I8, */
1246 : /* ***D 1 ' LX ='I6,' JMIN(LX) ='I8,' JMAX(LX) ='I8, */
1247 : /* ***D 2 ' LBIT(LX) ='I5,' NOV(LX) ='I8,' MISSLX(LX) =',I7) */
1248 :
1249 39979 : if (ktotal >= *nxy) {
1250 14 : goto L209;
1251 : }
1252 :
1253 : /* THE NEW GROUP A WILL BE THE PREVIOUS GROUP B. SET LIMITS, ETC. */
1254 :
1255 39965 : ibita = ibitb;
1256 39965 : mina = minb;
1257 39965 : maxa = maxb;
1258 39965 : minak = minbk;
1259 39965 : maxak = maxbk;
1260 39965 : mislla = misllb;
1261 39965 : nenda = nendb;
1262 39965 : kounta = kountb;
1263 39965 : ktotal += kounta;
1264 39965 : adda = FALSE_;
1265 39965 : goto L133;
1266 :
1267 : /* ************************************* */
1268 :
1269 : /* CALCULATE IBIT, THE NUMBER OF BITS NEEDED TO HOLD THE GROUP */
1270 : /* MINIMUM VALUES. */
1271 :
1272 : /* ************************************* */
1273 :
1274 14 : L209:
1275 14 : *ibit = 0;
1276 :
1277 14 : i__1 = *lx;
1278 39993 : for (l = 1; l <= i__1; ++l) {
1279 39979 : L210:
1280 40060 : if( *ibit == 31 )
1281 : {
1282 0 : *ier = -1;
1283 0 : goto L900;
1284 : }
1285 40060 : if (jmin[l] < ibxx2[*ibit]) {
1286 39979 : goto L220;
1287 : }
1288 81 : ++(*ibit);
1289 81 : goto L210;
1290 39979 : L220:
1291 : ;
1292 : }
1293 :
1294 : /* INSERT THE VALUE IN JMIN( ) TO BE USED FOR ALL MISSING */
1295 : /* VALUES WHEN LBIT( ) = 0. WHEN SECONDARY MISSING */
1296 : /* VALUES CAN BE PRESENT, LBIT(L) WILL NOT = 0. */
1297 :
1298 14 : if (*is523 == 1) {
1299 :
1300 4 : i__1 = *lx;
1301 39768 : for (l = 1; l <= i__1; ++l) {
1302 :
1303 39764 : if (lbit[l] == 0) {
1304 :
1305 3874 : if (misslx[l - 1] == *missp) {
1306 3874 : jmin[l] = ibxx2[*ibit] - 1;
1307 : }
1308 :
1309 : }
1310 :
1311 : /* L226: */
1312 : }
1313 :
1314 : }
1315 :
1316 : /* ************************************* */
1317 :
1318 : /* CALCULATE JBIT, THE NUMBER OF BITS NEEDED TO HOLD THE BITS */
1319 : /* NEEDED TO PACK THE VALUES IN THE GROUPS. BUT FIND AND */
1320 : /* REMOVE THE REFERENCE VALUE FIRST. */
1321 :
1322 : /* ************************************* */
1323 :
1324 : /* WRITE(KFILDO,228)CFEED,LX */
1325 : /* 228 FORMAT(A1,/' *****************************************' */
1326 : /* 1 /' THE GROUP WIDTHS LBIT( ) FOR ',I8,' GROUPS' */
1327 : /* 2 /' *****************************************') */
1328 : /* WRITE(KFILDO,229) (LBIT(J),J=1,MIN(LX,100)) */
1329 : /* 229 FORMAT(/' '20I6) */
1330 :
1331 14 : *lbitref = lbit[1];
1332 :
1333 14 : i__1 = *lx;
1334 39993 : for (k = 1; k <= i__1; ++k) {
1335 39979 : if (lbit[k] < *lbitref) {
1336 16 : *lbitref = lbit[k];
1337 : }
1338 : /* L230: */
1339 : }
1340 :
1341 14 : if (*lbitref != 0) {
1342 :
1343 3 : i__1 = *lx;
1344 28 : for (k = 1; k <= i__1; ++k) {
1345 25 : lbit[k] -= *lbitref;
1346 : /* L240: */
1347 : }
1348 :
1349 : }
1350 :
1351 : /* WRITE(KFILDO,241)CFEED,LBITREF */
1352 : /* 241 FORMAT(A1,/' *****************************************' */
1353 : /* 1 /' THE GROUP WIDTHS LBIT( ) AFTER REMOVING REFERENCE ', */
1354 : /* 2 I8, */
1355 : /* 3 /' *****************************************') */
1356 : /* WRITE(KFILDO,242) (LBIT(J),J=1,MIN(LX,100)) */
1357 : /* 242 FORMAT(/' '20I6) */
1358 :
1359 14 : *jbit = 0;
1360 :
1361 14 : i__1 = *lx;
1362 39993 : for (k = 1; k <= i__1; ++k) {
1363 39979 : L310:
1364 40015 : if (lbit[k] < ibxx2[*jbit]) {
1365 39979 : goto L320;
1366 : }
1367 36 : ++(*jbit);
1368 36 : goto L310;
1369 39979 : L320:
1370 : ;
1371 : }
1372 :
1373 : /* ************************************* */
1374 :
1375 : /* CALCULATE KBIT, THE NUMBER OF BITS NEEDED TO HOLD THE NUMBER */
1376 : /* OF VALUES IN THE GROUPS. BUT FIND AND REMOVE THE */
1377 : /* REFERENCE FIRST. */
1378 :
1379 : /* ************************************* */
1380 :
1381 : /* WRITE(KFILDO,321)CFEED,LX */
1382 : /* 321 FORMAT(A1,/' *****************************************' */
1383 : /* 1 /' THE GROUP SIZES NOV( ) FOR ',I8,' GROUPS' */
1384 : /* 2 /' *****************************************') */
1385 : /* WRITE(KFILDO,322) (NOV(J),J=1,MIN(LX,100)) */
1386 : /* 322 FORMAT(/' '20I6) */
1387 :
1388 14 : *novref = nov[1];
1389 :
1390 14 : i__1 = *lx;
1391 39993 : for (k = 1; k <= i__1; ++k) {
1392 39979 : if (nov[k] < *novref) {
1393 50 : *novref = nov[k];
1394 : }
1395 : /* L400: */
1396 : }
1397 :
1398 14 : if (*novref > 0) {
1399 :
1400 14 : i__1 = *lx;
1401 39993 : for (k = 1; k <= i__1; ++k) {
1402 39979 : nov[k] -= *novref;
1403 : /* L405: */
1404 : }
1405 :
1406 : }
1407 :
1408 : /* WRITE(KFILDO,406)CFEED,NOVREF */
1409 : /* 406 FORMAT(A1,/' *****************************************' */
1410 : /* 1 /' THE GROUP SIZES NOV( ) AFTER REMOVING REFERENCE ',I8, */
1411 : /* 2 /' *****************************************') */
1412 : /* WRITE(KFILDO,407) (NOV(J),J=1,MIN(LX,100)) */
1413 : /* 407 FORMAT(/' '20I6) */
1414 : /* WRITE(KFILDO,408)CFEED */
1415 : /* 408 FORMAT(A1,/' *****************************************' */
1416 : /* 1 /' THE GROUP REFERENCES JMIN( )' */
1417 : /* 2 /' *****************************************') */
1418 : /* WRITE(KFILDO,409) (JMIN(J),J=1,MIN(LX,100)) */
1419 : /* 409 FORMAT(/' '20I6) */
1420 :
1421 14 : *kbit = 0;
1422 :
1423 14 : i__1 = *lx;
1424 39993 : for (k = 1; k <= i__1; ++k) {
1425 39979 : L410:
1426 40077 : if (nov[k] < ibxx2[*kbit]) {
1427 39979 : goto L420;
1428 : }
1429 98 : ++(*kbit);
1430 98 : goto L410;
1431 39979 : L420:
1432 : ;
1433 : }
1434 :
1435 : /* DETERMINE WHETHER THE GROUP SIZES SHOULD BE REDUCED */
1436 : /* FOR SPACE EFFICIENCY. */
1437 :
1438 14 : if (ired == 0) {
1439 13 : reduce(kfildo, &jmin[1], &jmax[1], &lbit[1], &nov[1], lx, ndg, ibit,
1440 : jbit, kbit, novref, ibxx2, ier);
1441 :
1442 13 : if (*ier == 714 || *ier == 715) {
1443 : /* REDUCE HAS ABORTED. REEXECUTE PACK_GP WITHOUT REDUCE. */
1444 : /* PROVIDE FOR A NON FATAL RETURN FROM REDUCE. */
1445 1 : iersav = *ier;
1446 1 : ired = 1;
1447 1 : *ier = 0;
1448 1 : goto L102;
1449 : }
1450 :
1451 : }
1452 :
1453 13 : free(misslx);
1454 13 : misslx=0;
1455 :
1456 : /* CALL TIMPR(KFILDO,KFILDO,'END PACK_GP ') */
1457 13 : if (iersav != 0) {
1458 1 : *ier = iersav;
1459 1 : return 0;
1460 : }
1461 :
1462 : /* 900 IF(IER.NE.0)RETURN1 */
1463 :
1464 12 : L900:
1465 12 : free(misslx);
1466 12 : return 0;
1467 : } /* pack_gp__ */
1468 :
|
