1 #include <string.h>
2 #include <stdio.h>
3 #include <stdlib.h>
4 #include <time.h>
5 #include <sys/time.h>
6
7 /*
8 Every program which uses the CFITSIO interface must include the
9 the fitsio.h header file. This contains the prototypes for all
10 the routines and defines the error status values and other symbolic
11 constants used in the interface.
12 */
13 #include "fitsio.h"
14
15 #define minvalue(A,B) ((A) < (B) ? (A) : (B))
16
17 /* size of the image */
18 #define XSIZE 3000
19 #define YSIZE 3000
20
21 /* size of data buffer */
22 #define SHTSIZE 20000
23 static long sarray[ SHTSIZE ] = {SHTSIZE * 0};
24
25 /* no. of rows in binary table */
26 #define BROWS 2500000
27
28 /* no. of rows in ASCII table */
29 #define AROWS 400000
30
31 /* CLOCKS_PER_SEC should be defined by most compilers */
32 #if defined(CLOCKS_PER_SEC)
33 #define CLOCKTICKS CLOCKS_PER_SEC
34 #else
35 /* on SUN OS machine, CLOCKS_PER_SEC is not defined, so set its value */
36 #define CLOCKTICKS 1000000
37 #define difftime(A,B) ((double) A - (double) B)
38 #endif
39
40 /* define variables for measuring elapsed time */
41 clock_t scpu, ecpu;
42 time_t start, finish;
43 long startsec; /* start of elapsed time interval */
44 int startmilli; /* start of elapsed time interval */
45
46 int writeimage(fitsfile *fptr, int *status);
47 int writebintable(fitsfile *fptr, int *status);
48 int writeasctable(fitsfile *fptr, int *status);
49 int readimage(fitsfile *fptr, int *status);
50 int readatable(fitsfile *fptr, int *status);
51 int readbtable(fitsfile *fptr, int *status);
52 void printerror( int status);
53 int marktime(int *status);
54 int gettime(double *elapse, float *elapscpu, int *status);
55 int main(void);
56
main()57 int main()
58 {
59 /*************************************************************************
60 This program tests the speed of writing/reading FITS files with cfitsio
61 **************************************************************************/
62
63 FILE *diskfile;
64 fitsfile *fptr; /* pointer to the FITS file, defined in fitsio.h */
65 int status, ii;
66 long rawloop;
67 char filename[] = "speedcc.fit"; /* name for new FITS file */
68 char buffer[2880] = {2880 * 0};
69 time_t tbegin, tend;
70 float rate, size, elapcpu, cpufrac;
71 double elapse;
72
73 tbegin = time(0);
74
75 remove(filename); /* Delete old file if it already exists */
76
77 diskfile = fopen(filename,"w+b");
78 rawloop = XSIZE * YSIZE / 720;
79
80 printf(" ");
81 printf(" SIZE / ELAPSE(%%CPU) = RATE\n");
82 printf("RAW fwrite (2880 bytes/loop)... ");
83 marktime(&status);
84
85 for (ii = 0; ii < rawloop; ii++)
86 if (fwrite(buffer, 1, 2880, diskfile) != 2880)
87 printf("write error \n");
88
89 gettime(&elapse, &elapcpu, &status);
90
91 cpufrac = elapcpu / elapse * 100.;
92 size = 2880. * rawloop / 1000000.;
93 rate = size / elapse;
94 printf(" %4.1fMB/%6.3fs(%3.0f) = %5.2fMB/s\n", size, elapse, cpufrac,rate);
95
96 /* read back the binary records */
97 fseek(diskfile, 0, 0);
98
99 printf("RAW fread (2880 bytes/loop)... ");
100 marktime(&status);
101
102 for (ii = 0; ii < rawloop; ii++)
103 if (fread(buffer, 1, 2880, diskfile) != 2880)
104 printf("read error \n");
105
106 gettime(&elapse, &elapcpu, &status);
107
108 cpufrac = elapcpu / elapse * 100.;
109 size = 2880. * rawloop / 1000000.;
110 rate = size / elapse;
111 printf(" %4.1fMB/%6.3fs(%3.0f) = %5.2fMB/s\n", size, elapse, cpufrac,rate);
112
113 fclose(diskfile);
114 remove(filename);
115
116 status = 0;
117 fptr = 0;
118
119 if (fits_create_file(&fptr, filename, &status)) /* create new FITS file */
120 printerror( status);
121
122 if (writeimage(fptr, &status))
123 printerror( status);
124
125 if (writebintable(fptr, &status))
126 printerror( status);
127
128 if (writeasctable(fptr, &status))
129 printerror( status);
130
131 if (readimage(fptr, &status))
132 printerror( status);
133
134 if (readbtable(fptr, &status))
135 printerror( status);
136
137 if (readatable(fptr, &status))
138 printerror( status);
139
140 if (fits_close_file(fptr, &status))
141 printerror( status );
142
143 tend = time(0);
144 elapse = difftime(tend, tbegin) + 0.5;
145 printf("Total elapsed time = %.3fs, status = %d\n",elapse, status);
146 return(0);
147 }
148 /*--------------------------------------------------------------------------*/
writeimage(fitsfile * fptr,int * status)149 int writeimage(fitsfile *fptr, int *status)
150
151 /**************************************************/
152 /* write the primary array containing a 2-D image */
153 /**************************************************/
154 {
155 long nremain, ii;
156 float rate, size, elapcpu, cpufrac;
157 double elapse;
158
159 /* initialize FITS image parameters */
160 int bitpix = 32; /* 32-bit signed integer pixel values */
161 long naxis = 2; /* 2-dimensional image */
162 long naxes[2] = {XSIZE, YSIZE }; /* image size */
163
164 /* write the required keywords for the primary array image */
165 if ( fits_create_img(fptr, bitpix, naxis, naxes, status) )
166 printerror( *status );
167
168 printf("\nWrite %dx%d I*4 image, %d pixels/loop: ",XSIZE,YSIZE,SHTSIZE);
169 marktime(status);
170
171 nremain = XSIZE * YSIZE;
172 for (ii = 1; ii <= nremain; ii += SHTSIZE)
173 {
174 ffpprj(fptr, 0, ii, SHTSIZE, sarray, status);
175 }
176
177 ffflus(fptr, status); /* flush all buffers to disk */
178
179 gettime(&elapse, &elapcpu, status);
180
181 cpufrac = elapcpu / elapse * 100.;
182 size = XSIZE * 4. * YSIZE / 1000000.;
183 rate = size / elapse;
184 printf(" %4.1fMB/%6.3fs(%3.0f) = %5.2fMB/s\n", size, elapse, cpufrac,rate);
185
186 return( *status );
187 }
188 /*--------------------------------------------------------------------------*/
writebintable(fitsfile * fptr,int * status)189 int writebintable (fitsfile *fptr, int *status)
190
191 /*********************************************************/
192 /* Create a binary table extension containing 3 columns */
193 /*********************************************************/
194 {
195 int tfields = 2;
196 long nremain, ntodo, firstrow = 1, firstelem = 1, nrows;
197 float rate, size, elapcpu, cpufrac;
198 double elapse;
199
200 char extname[] = "Speed_Test"; /* extension name */
201
202 /* define the name, datatype, and physical units for the columns */
203 char *ttype[] = { "first", "second" };
204 char *tform[] = {"1J", "1J" };
205 char *tunit[] = { " ", " " };
206
207 /* append a new empty binary table onto the FITS file */
208
209 if ( fits_create_tbl( fptr, BINARY_TBL, BROWS, tfields, ttype, tform,
210 tunit, extname, status) )
211 printerror( *status );
212
213 /* get table row size and optimum number of rows to write per loop */
214 fits_get_rowsize(fptr, &nrows, status);
215 nrows = minvalue(nrows, SHTSIZE);
216 nremain = BROWS;
217
218 printf("Write %7drow x %dcol bintable %4ld rows/loop:", BROWS, tfields,
219 nrows);
220 marktime(status);
221
222 while(nremain)
223 {
224 ntodo = minvalue(nrows, nremain);
225 ffpclj(fptr, 1, firstrow, firstelem, ntodo, sarray, status);
226 ffpclj(fptr, 2, firstrow, firstelem, ntodo, sarray, status);
227 firstrow += ntodo;
228 nremain -= ntodo;
229 }
230
231 ffflus(fptr, status); /* flush all buffers to disk */
232
233 gettime(&elapse, &elapcpu, status);
234
235 cpufrac = elapcpu / elapse * 100.;
236 size = BROWS * 8. / 1000000.;
237 rate = size / elapse;
238 printf(" %4.1fMB/%6.3fs(%3.0f) = %5.2fMB/s\n", size, elapse, cpufrac,rate);
239
240 return( *status );
241 }
242 /*--------------------------------------------------------------------------*/
writeasctable(fitsfile * fptr,int * status)243 int writeasctable (fitsfile *fptr, int *status)
244
245 /*********************************************************/
246 /* Create an ASCII table extension containing 2 columns */
247 /*********************************************************/
248 {
249 int tfields = 2;
250 long nremain, ntodo, firstrow = 1, firstelem = 1;
251 long nrows;
252 float rate, size, elapcpu, cpufrac;
253 double elapse;
254
255 char extname[] = "Speed_Test"; /* extension name */
256
257 /* define the name, datatype, and physical units for the columns */
258 char *ttype[] = { "first", "second" };
259 char *tform[] = {"I6", "I6" };
260 char *tunit[] = { " ", " " };
261
262 /* append a new empty ASCII table onto the FITS file */
263 if ( fits_create_tbl( fptr, ASCII_TBL, AROWS, tfields, ttype, tform,
264 tunit, extname, status) )
265 printerror( *status );
266
267 /* get table row size and optimum number of rows to write per loop */
268 fits_get_rowsize(fptr, &nrows, status);
269 nrows = minvalue(nrows, SHTSIZE);
270 nremain = AROWS;
271
272 printf("Write %7drow x %dcol asctable %4ld rows/loop:", AROWS, tfields,
273 nrows);
274 marktime(status);
275
276 while(nremain)
277 {
278 ntodo = minvalue(nrows, nremain);
279 ffpclj(fptr, 1, firstrow, firstelem, ntodo, sarray, status);
280 ffpclj(fptr, 2, firstrow, firstelem, ntodo, sarray, status);
281 firstrow += ntodo;
282 nremain -= ntodo;
283 }
284
285 ffflus(fptr, status); /* flush all buffers to disk */
286
287 gettime(&elapse, &elapcpu, status);
288
289 cpufrac = elapcpu / elapse * 100.;
290 size = AROWS * 13. / 1000000.;
291 rate = size / elapse;
292 printf(" %4.1fMB/%6.3fs(%3.0f) = %5.2fMB/s\n", size, elapse, cpufrac,rate);
293
294 return( *status );
295 }
296 /*--------------------------------------------------------------------------*/
readimage(fitsfile * fptr,int * status)297 int readimage( fitsfile *fptr, int *status )
298
299 /*********************/
300 /* Read a FITS image */
301 /*********************/
302 {
303 int anynull, hdutype;
304 long nremain, ii;
305 long longnull = 0;
306 float rate, size, elapcpu, cpufrac;
307 double elapse;
308
309 /* move to the primary array */
310 if ( fits_movabs_hdu(fptr, 1, &hdutype, status) )
311 printerror( *status );
312
313 printf("\nRead back image ");
314 marktime(status);
315
316 nremain = XSIZE * YSIZE;
317 for (ii=1; ii <= nremain; ii += SHTSIZE)
318 {
319 ffgpvj(fptr, 0, ii, SHTSIZE, longnull, sarray, &anynull, status);
320 }
321
322 gettime(&elapse, &elapcpu, status);
323
324 cpufrac = elapcpu / elapse * 100.;
325 size = XSIZE * 4. * YSIZE / 1000000.;
326 rate = size / elapse;
327 printf(" %4.1fMB/%6.3fs(%3.0f) = %5.2fMB/s\n", size, elapse, cpufrac,rate);
328
329 return( *status );
330 }
331 /*--------------------------------------------------------------------------*/
readbtable(fitsfile * fptr,int * status)332 int readbtable( fitsfile *fptr, int *status )
333
334 /************************************************************/
335 /* read and print data values from the binary table */
336 /************************************************************/
337 {
338 int hdutype, anynull;
339 long nremain, ntodo, firstrow = 1, firstelem = 1;
340 long nrows;
341 long lnull = 0;
342 float rate, size, elapcpu, cpufrac;
343 double elapse;
344
345 /* move to the table */
346 if ( fits_movrel_hdu(fptr, 1, &hdutype, status) )
347 printerror( *status );
348
349 /* get table row size and optimum number of rows to read per loop */
350 fits_get_rowsize(fptr, &nrows, status);
351 nrows = minvalue(nrows, SHTSIZE);
352
353 /* read the columns */
354 nremain = BROWS;
355
356 printf("Read back BINTABLE ");
357 marktime(status);
358
359 while(nremain)
360 {
361 ntodo = minvalue(nrows, nremain);
362 ffgcvj(fptr, 1, firstrow, firstelem, ntodo,
363 lnull, sarray, &anynull, status);
364 ffgcvj(fptr, 2, firstrow, firstelem, ntodo,
365 lnull, sarray, &anynull, status);
366 firstrow += ntodo;
367 nremain -= ntodo;
368 }
369
370 gettime(&elapse, &elapcpu, status);
371
372 cpufrac = elapcpu / elapse * 100.;
373 size = BROWS * 8. / 1000000.;
374 rate = size / elapse;
375 printf(" %4.1fMB/%6.3fs(%3.0f) = %5.2fMB/s\n", size, elapse, cpufrac,rate);
376
377 return( *status );
378 }
379 /*--------------------------------------------------------------------------*/
readatable(fitsfile * fptr,int * status)380 int readatable( fitsfile *fptr, int *status )
381
382 /************************************************************/
383 /* read and print data values from an ASCII or binary table */
384 /************************************************************/
385 {
386 int hdutype, anynull;
387 long nremain, ntodo, firstrow = 1, firstelem = 1;
388 long nrows;
389 long lnull = 0;
390 float rate, size, elapcpu, cpufrac;
391 double elapse;
392
393 /* move to the table */
394 if ( fits_movrel_hdu(fptr, 1, &hdutype, status) )
395 printerror( *status );
396
397 /* get table row size and optimum number of rows to read per loop */
398 fits_get_rowsize(fptr, &nrows, status);
399 nrows = minvalue(nrows, SHTSIZE);
400
401 /* read the columns */
402 nremain = AROWS;
403
404 printf("Read back ASCII Table ");
405 marktime(status);
406
407 while(nremain)
408 {
409 ntodo = minvalue(nrows, nremain);
410 ffgcvj(fptr, 1, firstrow, firstelem, ntodo,
411 lnull, sarray, &anynull, status);
412 ffgcvj(fptr, 2, firstrow, firstelem, ntodo,
413 lnull, sarray, &anynull, status);
414 firstrow += ntodo;
415 nremain -= ntodo;
416 }
417
418 gettime(&elapse, &elapcpu, status);
419
420 cpufrac = elapcpu / elapse * 100.;
421 size = AROWS * 13. / 1000000.;
422 rate = size / elapse;
423 printf(" %4.1fMB/%6.3fs(%3.0f) = %5.2fMB/s\n", size, elapse, cpufrac,rate);
424
425 return( *status );
426 }
427 /*--------------------------------------------------------------------------*/
printerror(int status)428 void printerror( int status)
429 {
430 /*****************************************************/
431 /* Print out cfitsio error messages and exit program */
432 /*****************************************************/
433
434 char status_str[FLEN_STATUS], errmsg[FLEN_ERRMSG];
435
436 if (status)
437 fprintf(stderr, "\n*** Error occurred during program execution ***\n");
438
439 fits_get_errstatus(status, status_str); /* get the error description */
440 fprintf(stderr, "\nstatus = %d: %s\n", status, status_str);
441
442 /* get first message; null if stack is empty */
443 if ( fits_read_errmsg(errmsg) )
444 {
445 fprintf(stderr, "\nError message stack:\n");
446 fprintf(stderr, " %s\n", errmsg);
447
448 while ( fits_read_errmsg(errmsg) ) /* get remaining messages */
449 fprintf(stderr, " %s\n", errmsg);
450 }
451
452 exit( status ); /* terminate the program, returning error status */
453 }
454 /*--------------------------------------------------------------------------*/
marktime(int * status)455 int marktime( int *status)
456 {
457 double telapse;
458 time_t temp;
459 struct timeval tv;
460
461 temp = time(0);
462
463 /* Since elapsed time is only measured to the nearest second */
464 /* keep getting the time until the seconds tick just changes. */
465 /* This provides more consistent timing measurements since the */
466 /* intervals all start on an integer seconds. */
467
468 telapse = 0.;
469
470 scpu = clock();
471 start = time(0);
472 /*
473 while (telapse == 0.)
474 {
475 scpu = clock();
476 start = time(0);
477 telapse = difftime( start, temp );
478 }
479 */
480 gettimeofday (&tv, NULL);
481
482 startsec = tv.tv_sec;
483 startmilli = tv.tv_usec/1000;
484
485 return( *status );
486 }
487 /*--------------------------------------------------------------------------*/
gettime(double * elapse,float * elapscpu,int * status)488 int gettime(double *elapse, float *elapscpu, int *status)
489 {
490 struct timeval tv;
491 int stopmilli;
492 long stopsec;
493
494
495 gettimeofday (&tv, NULL);
496 ecpu = clock();
497 finish = time(0);
498
499 stopmilli = tv.tv_usec/1000;
500 stopsec = tv.tv_sec;
501
502
503 *elapse = (stopsec - startsec) + (stopmilli - startmilli)/1000.;
504
505 /* *elapse = difftime(finish, start) + 0.5; */
506 *elapscpu = (ecpu - scpu) * 1.0 / CLOCKTICKS;
507
508 return( *status );
509 }
510