From 8bb84e35d06fe6832ab92d7f0bd8b5e2a37745b3 Mon Sep 17 00:00:00 2001 From: radonnachie Date: Fri, 29 Jul 2022 09:48:12 +0000 Subject: [PATCH 1/9] *@ callback flag_file_output, enum *_file_format_t --- rawspec.c | 68 ++++++++++++++++++++------------ rawspec_callback.h | 11 +++++- rawspec_file.c | 97 ++++++++++++++++++++++++++++------------------ 3 files changed, 112 insertions(+), 64 deletions(-) diff --git a/rawspec.c b/rawspec.c index 73d9084..612a214 100644 --- a/rawspec.c +++ b/rawspec.c @@ -213,7 +213,7 @@ int main(int argc, char *argv[]) int flag_debugging = 0; // FBH5 fields - int flag_fbh5_output = 0; + enum rawspec_callback_file_format_t flag_file_output = FILE_FORMAT_FBSIGPROC; // For net data rate rate calculations double rate = 6.0; @@ -245,7 +245,7 @@ int main(int argc, char *argv[]) break; case 'j': // FBH5 output format requested - flag_fbh5_output = 1; + flag_file_output = FILE_FORMAT_FBH5; break; case 'z': // Selected dynamic debugging @@ -403,10 +403,14 @@ int main(int argc, char *argv[]) // If writing output files, show the format used if(output_mode == RAWSPEC_FILE) { - if(flag_fbh5_output) - printf("writing output files in FBH5 format\n"); - else - printf("writing output files in SIGPROC Filterbank format\n"); + switch(flag_file_output) { + case FILE_FORMAT_FBH5: + printf("writing output files in FBH5 format\n"); + break; + case FILE_FORMAT_FBSIGPROC: + printf("writing output files in SIGPROC Filterbank format\n"); + break; + } } // If schan is non-zero, nchan must be too @@ -495,12 +499,13 @@ int main(int argc, char *argv[]) // Init callback file descriptors to sentinal values cb_data[i].fd = malloc(sizeof(int)); cb_data[i].fd[0] = -1; - if(flag_fbh5_output) { - cb_data[i].flag_fbh5_output = 1; + + cb_data[i].flag_file_output = flag_file_output; + switch(flag_file_output) { + case FILE_FORMAT_FBH5: cb_data[i].fbh5_ctx_ant = malloc(sizeof(fbh5_context_t)); cb_data[i].fbh5_ctx_ant[0].active = 0; - } else { - cb_data[i].flag_fbh5_output = 0; + break; } } @@ -639,12 +644,14 @@ int main(int argc, char *argv[]) // For each antenna ..... for(j=0; jflag_fbh5_output) + switch(cb_data->flag_file_output) { + case FILE_FORMAT_FBH5: strcpy(fileext, "h5"); - else + break; + case FILE_FORMAT_FBSIGPROC: strcpy(fileext, "fil"); + break; + } if(dest && dest[0]) { // Look for last '/' in stem basename = strrchr(stem, '/'); @@ -40,7 +44,9 @@ int open_output_file(callback_data_t *cb_data, const char * dest, const char *st snprintf(fname, PATH_MAX, "%s.rawspec.%04d.%s", stem, output_idx, fileext); } fname[PATH_MAX] = '\0'; - if(cb_data->flag_fbh5_output) { + + switch(cb_data->flag_file_output) { + case FILE_FORMAT_FBH5: // Open an FBH5 output file. // If antenna_index < 0, then use the ICS context; // Else, use the indicated antenna context. @@ -55,17 +61,19 @@ int open_output_file(callback_data_t *cb_data, const char * dest, const char *st if(cb_data->debug_callback) printf("open_output_file: fbh5_open(%s) successful\n", fname); return ENABLER_FD_FOR_FBH5; - } - // Open a SIGPROC Filterbank output file. - fd = open(fname, O_WRONLY | O_CREAT | O_TRUNC, 0664); - if(fd == -1) { - cb_data->exit_soon = 1; - perror(fname); - } else { - posix_fadvise(fd, 0, 0, POSIX_FADV_DONTNEED); + case FILE_FORMAT_FBSIGPROC: + // Open a SIGPROC Filterbank output file. + fd = open(fname, O_WRONLY | O_CREAT | O_TRUNC, 0664); + if(fd == -1) { + cb_data->exit_soon = 1; + perror(fname); + } else { + posix_fadvise(fd, 0, 0, POSIX_FADV_DONTNEED); + } + return fd; + } - return fd; } // Open one or more output Filterbank files for the following cases: @@ -93,8 +101,11 @@ int open_output_file_per_antenna_and_write_header(callback_data_t *cb_data, cons } // Write filterbank header to output file if SIGPROC. - if(! cb_data->flag_fbh5_output) + switch(cb_data->flag_file_output) { + case FILE_FORMAT_FBSIGPROC: fb_fd_write_header(cb_data->fd[i], &cb_data->fb_hdr); + break; + } } return 0; } @@ -108,25 +119,29 @@ void * dump_file_thread_func(void *arg) if(cb_data->fd && cb_data->h_pwrbuf && (cb_data->Nant == 1)) { if(cb_data->debug_callback) printf("dump_file_thread_func: write for nants=0\n"); - if(cb_data->flag_fbh5_output) { - retcode = fbh5_write(&(cb_data->fbh5_ctx_ant[0]), - &(cb_data->fb_hdr), - cb_data->h_pwrbuf, - cb_data->h_pwrbuf_size, - cb_data->debug_callback); - if(retcode != 0) { - cb_data->exit_soon = 1; - cb_data->output_thread_valid = 0; - } - } else { // SIGPROC Filterbank - if(write(cb_data->fd[0], - cb_data->h_pwrbuf, - cb_data->h_pwrbuf_size) < 0) { - cb_data->exit_soon = 1; - cb_data->output_thread_valid = 0; - fprintf(stderr, "SIGPROC-WRITE-ERROR\n"); - } - } // if(cb_data->flag_fbh5_output) + + switch(cb_data->flag_file_output) { + case FILE_FORMAT_FBH5: + retcode = fbh5_write(&(cb_data->fbh5_ctx_ant[0]), + &(cb_data->fb_hdr), + cb_data->h_pwrbuf, + cb_data->h_pwrbuf_size, + cb_data->debug_callback); + if(retcode != 0) { + cb_data->exit_soon = 1; + cb_data->output_thread_valid = 0; + } + break; + case FILE_FORMAT_FBSIGPROC: + if(write(cb_data->fd[0], + cb_data->h_pwrbuf, + cb_data->h_pwrbuf_size) < 0) { + cb_data->exit_soon = 1; + cb_data->output_thread_valid = 0; + fprintf(stderr, "SIGPROC-WRITE-ERROR\n"); + } + break; + } } // Multiple antennas, split output @@ -145,7 +160,9 @@ void * dump_file_thread_func(void *arg) } if(cb_data->debug_callback) printf("dump_file_thread_func: write for antenna %ld-%ld-%ld\n", k, j, i); - if(cb_data->flag_fbh5_output) { + + switch(cb_data->flag_file_output) { + case FILE_FORMAT_FBH5: retcode = fbh5_write(&(cb_data->fbh5_ctx_ant[i]), &(cb_data->fb_hdr), cb_data->h_pwrbuf + i * ant_stride + j * pol_stride + k * spectra_stride, @@ -155,7 +172,8 @@ void * dump_file_thread_func(void *arg) cb_data->exit_soon = 1; cb_data->output_thread_valid = 0; } - } else { // SIGPROC Filterbank + break; + case FILE_FORMAT_FBSIGPROC: if(write(cb_data->fd[i], cb_data->h_pwrbuf + i * ant_stride + j * pol_stride + k * spectra_stride, ant_stride * sizeof(float)) < 0) { @@ -163,7 +181,8 @@ void * dump_file_thread_func(void *arg) cb_data->output_thread_valid = 0; fprintf(stderr, "SIGPROC-WRITE-ERROR\n"); } - } // if(cb_data->flag_fbh5_output) + break; + } } // for(size_t i = 0; i < cb_data->Nant; i++) } // for(size_t j = 0; j < cb_data->fb_hdr.nifs; j++) } // for(size_t k = 0; k < cb_data->Nds; k++) @@ -174,7 +193,9 @@ void * dump_file_thread_func(void *arg) if(cb_data->fd_ics && cb_data->h_icsbuf) { if(cb_data->debug_callback) printf("dump_file_thread_func: write for ICS\n"); - if(cb_data->flag_fbh5_output) { + + switch(cb_data->flag_file_output) { + case FILE_FORMAT_FBH5: retcode = fbh5_write(&(cb_data->fbh5_ctx_ics), &(cb_data->fb_hdr), cb_data->h_icsbuf, @@ -184,7 +205,8 @@ void * dump_file_thread_func(void *arg) cb_data->exit_soon = 1; cb_data->output_thread_valid = 0; } - } else { // SIGPROC Filterbank + break; + case FILE_FORMAT_FBSIGPROC: if(write(cb_data->fd_ics, cb_data->h_icsbuf, cb_data->h_pwrbuf_size/cb_data->Nant) < 0) { @@ -192,6 +214,7 @@ void * dump_file_thread_func(void *arg) cb_data->output_thread_valid = 0; fprintf(stderr, "SIGPROC-WRITE-ERROR\n"); } + break; } } From 0daf3fd25c93ab61c50b1b938ed7e20b7761ce9d Mon Sep 17 00:00:00 2001 From: radonnachie Date: Fri, 29 Jul 2022 15:09:37 +0000 Subject: [PATCH 2/9] *= rawspec_raw_parse_header, guppirawc99 dependency --- Makefile | 41 +- hget.h | 2544 -------------------------------------------- rawspec_rawutils.c | 110 +- rawspec_rawutils.h | 2 + 4 files changed, 119 insertions(+), 2578 deletions(-) delete mode 100644 hget.h diff --git a/Makefile b/Makefile index 9cbae24..b8e3e79 100644 --- a/Makefile +++ b/Makefile @@ -39,6 +39,39 @@ LINKH5:= $(LIBDIR_H5) -l $(LIBHDF5) -l $(LIBHDF5_HL) # End HDF5 definitions +# Begin GuppiRawC99 definitions + +# if we have pkg-config installed, do we have guppirawc99 libraries installed? +ifeq ($(HAVEPKG),yes) + HAVEGRC99 = $(shell pkg-config --exists guppirawc99 && echo yes) +endif + +# depending on the above questions use, in order of priority: +# 1. User supplied values for INCDIR_GRC99 and LIBDIR_GRC99 +# 2. Whatever pkg-config supplies for those variables +# 3. Failing either the above, default values /usr/local/include and /usr/local/lib + +ifeq ($(INCDIR_GRC99),) + ifeq ($(HAVEGRC99),yes) + INCDIR_GRC99 = $(shell pkg-config --cflags guppirawc99) + else + INCDIR_GRC99 = -I/usr/local/include + endif +endif + +ifeq ($(LIBDIR_GRC99),) + ifeq ($(HAVEGRC99),yes) + LIBDIR_GRC99 = $(shell pkg-config --libs guppirawc99) + else + LIBDIR_GRC99 = -L/usr/local/lib + endif +endif + +LIBGRC99= :libguppirawc99.so +LINKGRC99:= $(LIBDIR_GRC99) -l $(LIBGRC99) + +# End GuppiRawC99 definitions + CUDA_DIR ?= $(CUDA_ROOT) CUDA_PATH ?= $(CUDA_DIR) @@ -47,7 +80,7 @@ CXX ?= g++ HOST_COMPILER ?= $(CXX) NVCC := $(CUDA_PATH)/bin/nvcc -ccbin $(HOST_COMPILER) -CFLAGS = -ggdb -fPIC -I$(CUDA_PATH)/include $(INCDIR_H5) +CFLAGS = -ggdb -fPIC -I$(CUDA_PATH)/include $(INCDIR_H5) $(INCDIR_GRC99) ifdef DEBUG_CALLBACKS CFLAGS += -DDEBUG_CALLBACKS=$(DEBUG_CALLBACKS) endif @@ -98,7 +131,7 @@ rawspec_gpu.o: rawspec.h rawspec_version.h cufft_error_name.h rawspec_socket.o: rawspec_socket.h rawspec.h \ rawspec_callback.h rawspec_fbutils.h rawspectest.o: rawspec.h -rawspec_rawutils.o: rawspec_rawutils.h hget.h +rawspec_rawutils.o: rawspec_rawutils.h # Begin fbh5 objects fbh5_open.o: fbh5_defs.h rawspec_callback.h rawspec_fbutils.h @@ -111,11 +144,11 @@ fbh5_util.o: fbh5_defs.h rawspec_callback.h rawspec_fbutils.h $(VERBOSE) $(NVCC) $(NVCC_FLAGS) -dc $(GENCODE_FLAGS) -o $@ -c $< librawspec.so: rawspec_gpu.o rawspec_fbutils.o rawspec_rawutils.o fbh5_open.o fbh5_close.o fbh5_write.o fbh5_util.o - $(VERBOSE) $(NVCC) -shared $(NVCC_FLAGS) $(GENCODE_FLAGS) -o $@ $^ $(CUDA_STATIC_LIBS) $(LINKH5) + $(VERBOSE) $(NVCC) -shared $(NVCC_FLAGS) $(GENCODE_FLAGS) -o $@ $^ $(CUDA_STATIC_LIBS) $(LINKH5) $(LINKGRC99) rawspec: librawspec.so rawspec: rawspec.o rawspec_file.o rawspec_socket.o - $(VERBOSE) $(NVCC) $(NVCC_FLAGS) $(GENCODE_FLAGS) -o $@ $^ -L. -lrawspec $(LINKH5) + $(VERBOSE) $(NVCC) $(NVCC_FLAGS) $(GENCODE_FLAGS) -o $@ $^ -L. -lrawspec $(LINKH5) $(LINKGRC99) rawspectest: librawspec.so rawspectest: rawspectest.o diff --git a/hget.h b/hget.h deleted file mode 100644 index 1d40ef0..0000000 --- a/hget.h +++ /dev/null @@ -1,2544 +0,0 @@ -// This file was originally hget.c by Doug Mink of the SAO. It has been -// converted into a .h file of static functions for internal use within a C -// compilation unit. You may want to consider using the original hget.c/hput.c -// rather than this file. - -#ifndef HGET_H -#define HGET_H - -#include -#include -#include -#include /* NULL, strlen, strstr, strcpy */ -#include - -#define VLENGTH 81 -#define multiline 0 -#define lhead0 0 - -static char * -hgetc(const char *hstring, const char *keyword0, char *value_buffer); - -static int -hgets(const char *hstring, const char *keyword, const const int lstr, char *str); - -static int -isnum(const char *string); - -static char * -ksearch(const char *hstring, const char *keyword); - -static char * -strsrch(const char *s1, const char *s2); - -static char * -strnsrch(const char *s1, const char *s2, const const int ls1); - -static char * -strcsrch(const char *s1, const char *s2); - -static char * -strncsrch(const char *s1, const char *s2, const const int ls1); - -/*** File libwcs/hget.c - *** August 22, 2007 - *** By Doug Mink, dmink@cfa.harvard.edu - *** Harvard-Smithsonian Center for Astrophysics - *** Copyright (C) 1994-2007 - *** Smithsonian Astrophysical Observatory, Cambridge, MA, USA - - This library is free software; you can redistribute it and/or - modify it under the terms of the GNU Lesser General Public - License as published by the Free Software Foundation; either - version 2 of the License, or (at your option) any later version. - - This library is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU - Lesser General Public License for more details. - - You should have received a copy of the GNU Lesser General Public - License along with this library; if not, write to the Free Software - Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA - - Correspondence concerning WCSTools should be addressed as follows: - Internet email: dmink@cfa.harvard.edu - Postal address: Doug Mink - Smithsonian Astrophysical Observatory - 60 Garden St. - Cambridge, MA 02138 USA - - * Module: hget.c (Get FITS Header parameter values) - * Purpose: Extract values for variables from FITS header string - * Subroutine: hgeti2 (hstring,keyword,ival) returns short integer - * Subroutine: hgeti4c (hstring,keyword,wchar,ival) returns long integer - * Subroutine: hgeti4 (hstring,keyword,ival) returns long integer - * Subroutine: hgetr4 (hstring,keyword,rval) returns real - * Subroutine: hgetra (hstring,keyword,ra) returns double RA in degrees - * Subroutine: hgetdec (hstring,keyword,dec) returns double Dec in degrees - * Subroutine: hgetr8c (hstring,keyword,wchar,dval) returns double - * Subroutine: hgetr8 (hstring,keyword,dval) returns double - * Subroutine: hgetl (hstring,keyword,lval) returns logical int (0=F, 1=T) - * Subroutine: hgetsc (hstring,keyword,wchar,lstr,str) returns character string - * Subroutine: hgets (hstring,keyword, lstr, str) returns character string - * Subroutine: hgetm (hstring,keyword, lstr, str) returns multi-keyword string - * Subroutine: hgetdate (hstring,keyword,date) returns date as fractional year - * Subroutine: hgetndec (hstring, keyword, ndec) returns number of dec. places - * Subroutine: hgetc (hstring,keyword,value_buffer) returns character string - * Subroutine: blsearch (hstring,keyword) returns pointer to blank lines - before keyword - * Subroutine: ksearch (hstring,keyword) returns pointer to header string entry - * Subroutine: str2ra (in) converts string to right ascension in degrees - * Subroutine: str2dec (in) converts string to declination in degrees - * Subroutine: strsrch (s1, s2) finds string s2 in null-terminated string s1 - * Subroutine: strnsrch (s1, s2, ls1) finds string s2 in ls1-byte string s1 - * Subroutine: hlength (header,lhead) sets length of FITS header for searching - * Subroutine: isnum (string) returns 1 if integer, 2 if fp number, else 0 - * Subroutine: notnum (string) returns 0 if number, else 1 - * Subroutine: numdec (string) returns number of decimal places in numeric string - * Subroutine: strfix (string,blankfill,zerodrop) removes extraneous characters - */ - -#if 0 -#include /* NULL, strlen, strstr, strcpy */ -#include -#include "fitshead.h" /* FITS header extraction subroutines */ -#include - -#ifndef VMS - -#include - -#else -#define INT_MAX 2147483647 /* Biggest number that can fit in long */ -#define SHRT_MAX 32767 -#endif -#define VLENGTH 81 - -static char *hgetc(const char *hstring, const char *keyword0, char *value_buffer); - -/* Don't know why this is global, seems to be just a - * temp buffer for parsing strings. Made local copies in - * each function that uses it. --PBD - */ -//static char val[VLENGTH+1]; - -/* See comments at start of hgetm. Changed multiline to const 0. - * This is fine if we never need multiline strings. - * --PBD - */ -//static int multiline = 0; -const static int multiline = 0; - -/* This looks like a hack to get around header strings that - * aren't null terminated. We can fix it at zero as long as - * we always have a null somewhere before our header memory runs - * out. --PBD - */ -//static int lhead0 = 0; /* Length of header string */ -const static int lhead0 = 0; /* Length of header string */ - -/* Set the length of the header string, if not terminated by NULL */ -int -hlength(const char *header, int lhead) -/* FITS header */ -/* Maximum length of FITS header */ -{ - char *hend; -#if 0 - if (lhead0 > 0) - lhead0 = lhead; - else { - lhead0 = 0; - hend = ksearch (header,"END"); - lhead0 = hend + 80 - header; - } -#endif - /* Always return header length based on END. --PBD*/ - hend = ksearch(header, "END"); - return (hend + 80 - header); -} - -/* Return the length of the header string, computing it if lhead0 not set */ -int -gethlength(char *header) -/* FITS header */ -{ - if (lhead0 > 0) - { - return (lhead0); - } - else - { - return (hlength(header, 0)); - } -} - -/* Extract Integer*4 value for variable from FITS header string */ - -int -hgeti4c(const char *hstring, const char *keyword, const char *wchar, int *ival) - -/* character string containing FITS header information - in the format = {/ } */ -/* character string containing the name of the keyword - the value of which is returned. hget searches for - a line beginning with this string. if "[n]" is - present, the n'th token in the value is returned. - (the first 8 characters must be unique) */ -/* Character of multiple WCS header; =0 if unused */ -/* Keyword value returned */ -{ - char keyword1[16]; - int lkey; - - if (wchar[0] < (char) 64) - { - return (hgeti4(hstring, keyword, ival)); - } - else - { - strcpy(keyword1, keyword); - lkey = strlen(keyword); - keyword1[lkey] = wchar[0]; - keyword1[lkey + 1] = (char) 0; - return (hgeti4(hstring, keyword1, ival)); - } -} -#endif // 0 - -/* Extract integer*8 value for variable from FITS header string */ - -static int -hgeti8 (hstring,keyword,i8val) - -const char *hstring; /* character string containing FITS header information - in the format = {/ } */ -const char *keyword; /* character string containing the name of the keyword - the value of which is returned. hget searches for a - line beginning with this string. if "[n]" is present, - the n'th token in the value is returned. - (the first 8 characters must be unique) */ -int64_t *i8val; -{ - char *value; - char *endptr; - char value_buffer[VLENGTH + 1]; - - /* Get value and comment from header string */ - value = hgetc(hstring, keyword, value_buffer); - - /* Translate value from ASCII to binary */ - if (value != NULL) { - if (value[0] == '#') value++; - *i8val = strtoll(value, &endptr, 0); - if(endptr && endptr[0]) { - fprintf(stderr, "%s:%s got invalid integer character '%c' (%d)\n", - __FUNCTION__, keyword, endptr[0], endptr[0]); - *i8val = (long long)atof(value); - } - return 1; - } else { - return 0; - } -} - -/* Extract unsigned integer*8 value for variable from FITS header string */ - -static int -hgetu8 (hstring,keyword,i8val) - -const char *hstring; /* character string containing FITS header information - in the format = {/ } */ -const char *keyword; /* character string containing the name of the keyword - the value of which is returned. hget searches for a - line beginning with this string. if "[n]" is present, - the n'th token in the value is returned. - (the first 8 characters must be unique) */ -uint64_t *i8val; -{ - char *value; - char *endptr; - char value_buffer[VLENGTH + 1]; - - /* Get value and comment from header string */ - value = hgetc (hstring,keyword, value_buffer); - - /* Translate value from ASCII to binary */ - if (value != NULL) { - if (value[0] == '#') value++; - *i8val = strtoull(value, &endptr, 0); - if(endptr && endptr[0]) { - fprintf(stderr, "%s:%s got invalid integer character '%c' (%d)\n", - __FUNCTION__, keyword, endptr[0], endptr[0]); - *i8val = (unsigned long long)atof(value); - } - return 1; - } else { - return 0; - } -} - -/* Extract long value for variable from FITS header string */ - -static int -hgeti4(const char *hstring, const char *keyword, int *ival) - -/* character string containing FITS header information -in the format = {/ } */ -/* character string containing the name of the keyword -the value of which is returned. hget searches for a -line beginning with this string. if "[n]" is present, -the n'th token in the value is returned. -(the first 8 characters must be unique) */ - -{ - char *value; - double dval; - int minint; - int lval; - char *dchar; - char val[VLENGTH + 1]; - char value_buffer[VLENGTH + 1]; - - /* Get value and comment from header string */ - value = hgetc(hstring, keyword, value_buffer); - - /* Translate value from ASCII to binary */ - if (value != NULL) - { - if (value[0] == '#') - { - value++; - } - minint = -INT_MAX - 1; - lval = strlen(value); - if (lval > VLENGTH) - { - strncpy(val, value, VLENGTH); - val[VLENGTH] = (char) 0; - } - else - { - strcpy(val, value); - } - if (isnum(val) == 2) - { - if ((dchar = strchr(val, 'D'))) - { - *dchar = 'e'; - } - if ((dchar = strchr(val, 'd'))) - { - *dchar = 'e'; - } - if ((dchar = strchr(val, 'E'))) - { - *dchar = 'e'; - } - } - dval = atof(val); - if (dval + 0.001 > INT_MAX) - { - *ival = INT_MAX; - } - else if (dval >= 0) - { - *ival = (int) (dval + 0.001); - } - else if (dval - 0.001 < minint) - { - *ival = minint; - } - else - { - *ival = (int) (dval - 0.001); - } - return (1); - } - else - { - return (0); - } -} - -/* Extract unsigned integer*4 value for variable from FITS header string */ - -static uint32_t -hgetu4 (hstring,keyword,i4val) - -const char *hstring; /* character string containing FITS header information - in the format = {/ } */ -const char *keyword; /* character string containing the name of the keyword - the value of which is returned. hget searches for a - line beginning with this string. if "[n]" is present, - the n'th token in the value is returned. - (the first 8 characters must be unique) */ -uint32_t *i4val; -{ - char *value; - char *endptr; - char value_buffer[VLENGTH + 1]; - - /* Get value and comment from header string */ - value = hgetc (hstring,keyword, value_buffer); - - /* Translate value from ASCII to binary */ - if (value != NULL) { - if (value[0] == '#') value++; - *i4val = strtoul(value, &endptr, 0); - if(endptr && endptr[0]) { - fprintf(stderr, "%s:%s got invalid integer character '%c' (%d)\n", - __FUNCTION__, keyword, endptr[0], endptr[0]); - *i4val = (uint32_t)atof(value); - } - return 1; - } else { - return 0; - } -} -#if 0 - -/* Extract integer*2 value for variable from fits header string */ - -int -hgeti2(const char *hstring, const char *keyword, short int *ival) - -/* character string containing FITS header information -in the format = {/ } */ -/* character string containing the name of the keyword -the value of which is returned. hget searches for a -line beginning with this string. if "[n]" is present, -the n'th token in the value is returned. -(the first 8 characters must be unique) */ - -{ - char *value; - double dval; - int minshort; - int lval; - char *dchar; - char val[VLENGTH + 1]; - char value_buffer[VLENGTH + 1]; - - /* Get value and comment from header string */ - value = hgetc(hstring, keyword, value_buffer); - - /* Translate value from ASCII to binary */ - if (value != NULL) - { - if (value[0] == '#') - { - value++; - } - lval = strlen(value); - if (lval > VLENGTH) - { - strncpy(val, value, VLENGTH); - val[VLENGTH] = (char) 0; - } - else - { - strcpy(val, value); - } - if (isnum(val) == 2) - { - if ((dchar = strchr(val, 'D'))) - { - *dchar = 'e'; - } - if ((dchar = strchr(val, 'd'))) - { - *dchar = 'e'; - } - if ((dchar = strchr(val, 'E'))) - { - *dchar = 'e'; - } - } - dval = atof(val); - minshort = -SHRT_MAX - 1; - if (dval + 0.001 > SHRT_MAX) - { - *ival = SHRT_MAX; - } - else if (dval >= 0) - { - *ival = (short) (dval + 0.001); - } - else if (dval - 0.001 < minshort) - { - *ival = minshort; - } - else - { - *ival = (short) (dval - 0.001); - } - return (1); - } - else - { - return (0); - } -} - -/* Extract real value for variable from FITS header string */ - -int -hgetr4(const char *hstring, const char *keyword, float *rval) - -/* character string containing FITS header information -in the format = {/ } */ -/* character string containing the name of the keyword -the value of which is returned. hget searches for a -line beginning with this string. if "[n]" is present, -the n'th token in the value is returned. -(the first 8 characters must be unique) */ - -{ - char *value; - int lval; - char *dchar; - char val[VLENGTH + 1]; - char value_buffer[VLENGTH + 1]; - - /* Get value and comment from header string */ - value = hgetc(hstring, keyword, value_buffer); - - /* translate value from ASCII to binary */ - if (value != NULL) - { - if (value[0] == '#') - { - value++; - } - lval = strlen(value); - if (lval > VLENGTH) - { - strncpy(val, value, VLENGTH); - val[VLENGTH] = (char) 0; - } - else - { - strcpy(val, value); - } - if (isnum(val) == 2) - { - if ((dchar = strchr(val, 'D'))) - { - *dchar = 'e'; - } - if ((dchar = strchr(val, 'd'))) - { - *dchar = 'e'; - } - if ((dchar = strchr(val, 'E'))) - { - *dchar = 'e'; - } - } - *rval = (float) atof(val); - return (1); - } - else - { - return (0); - } -} - - -/* Extract real*8 right ascension in degrees from FITS header string */ - -int -hgetra(const char *hstring, const char *keyword, double *dval) - -/* character string containing FITS header information -in the format = {/ } */ -/* character string containing the name of the keyword -the value of which is returned. hget searches for a -line beginning with this string. if "[n]" is present, -the n'th token in the value is returned. -(the first 8 characters must be unique) */ -/* Right ascension in degrees (returned) */ -{ - char *value; - char value_buffer[VLENGTH + 1]; - - /* Get value from header string */ - value = hgetc(hstring, keyword, value_buffer); - - /* Translate value from ASCII colon-delimited string to binary */ - if (value != NULL) - { - *dval = str2ra(value); - return (1); - } - else - { - return (0); - } -} - - -/* Extract real*8 declination in degrees from FITS header string */ - -int -hgetdec(const char *hstring, const char *keyword, double *dval) - -/* character string containing FITS header information -in the format = {/ } */ -/* character string containing the name of the keyword -the value of which is returned. hget searches for a -line beginning with this string. if "[n]" is present, -the n'th token in the value is returned. -(the first 8 characters must be unique) */ -/* Right ascension in degrees (returned) */ -{ - char *value; - char value_buffer[VLENGTH + 1]; - - /* Get value from header string */ - value = hgetc(hstring, keyword, value_buffer); - - /* Translate value from ASCII colon-delimited string to binary */ - if (value != NULL) - { - *dval = str2dec(value); - return (1); - } - else - { - return (0); - } -} - - -/* Extract real*8 value for variable from FITS header string */ - -int -hgetr8c(const char *hstring, const char *keyword, const char *wchar, double *dval) - -/* character string containing FITS header information - in the format = {/ } */ -/* character string containing the name of the keyword - the value of which is returned. hget searches for - a line beginning with this string. if "[n]" is - present, the n'th token in the value is returned. - (the first 8 characters must be unique) */ -/* Character of multiple WCS header; =0 if unused */ -/* Keyword value returned */ -{ - char keyword1[16]; - int lkey; - - if (wchar[0] < (char) 64) - { - return (hgetr8(hstring, keyword, dval)); - } - else - { - strcpy(keyword1, keyword); - lkey = strlen(keyword); - keyword1[lkey] = wchar[0]; - keyword1[lkey + 1] = (char) 0; - return (hgetr8(hstring, keyword1, dval)); - } -} -#endif // 0 - - -/* Extract real*8 value for variable from FITS header string */ - -static int -hgetr8(const char *hstring, const char *keyword, double *dval) - -/* character string containing FITS header information -in the format = {/ } */ -/* character string containing the name of the keyword -the value of which is returned. hget searches for a -line beginning with this string. if "[n]" is present, -the n'th token in the value is returned. -(the first 8 characters must be unique) */ - -{ - char *value; - int lval; - char *dchar; - char val[VLENGTH + 1]; - char value_buffer[VLENGTH + 1]; - - /* Get value and comment from header string */ - value = hgetc(hstring, keyword, value_buffer); - - /* Translate value from ASCII to binary */ - if (value != NULL) - { - if (value[0] == '#') - { - value++; - } - lval = strlen(value); - if (lval > VLENGTH) - { - strncpy(val, value, VLENGTH); - val[VLENGTH] = (char) 0; - } - else - { - strcpy(val, value); - } - if (isnum(val) == 2) - { - if ((dchar = strchr(val, 'D'))) - { - *dchar = 'e'; - } - if ((dchar = strchr(val, 'd'))) - { - *dchar = 'e'; - } - if ((dchar = strchr(val, 'E'))) - { - *dchar = 'e'; - } - } - *dval = atof(val); - return (1); - } - else - { - return (0); - } -} -#if 0 - -/* Extract logical value for variable from FITS header string */ - -int -hgetl(const char *hstring, const char *keyword, int *ival) - -/* character string containing FITS header information -in the format = {/ } */ -/* character string containing the name of the keyword -the value of which is returned. hget searches for a -line beginning with this string. if "[n]" is present, -the n'th token in the value is returned. -(the first 8 characters must be unique) */ - -{ - char *value; - char newval; - int lval; - char val[VLENGTH + 1]; - char value_buffer[VLENGTH + 1]; - - /* Get value and comment from header string */ - value = hgetc(hstring, keyword, value_buffer); - - /* Translate value from ASCII to binary */ - if (value != NULL) - { - lval = strlen(value); - if (lval > VLENGTH) - { - strncpy(val, value, VLENGTH); - val[VLENGTH] = (char) 0; - } - else - { - strcpy(val, value); - } - newval = val[0]; - if (newval == 't' || newval == 'T') - { - *ival = 1; - } - else - { - *ival = 0; - } - return (1); - } - else - { - return (0); - } -} - - -/* Extract real*8 date from FITS header string (dd/mm/yy or dd-mm-yy) */ - -int -hgetdate(const char *hstring, const char *keyword, double *dval) - -/* character string containing FITS header information -in the format = {/ } */ -/* character string containing the name of the keyword -the value of which is returned. hget searches for a -line beginning with this string. if "[n]" is present, -the n'th token in the value is returned. -(the first 8 characters must be unique) */ - -{ - double yeardays, seconds, fday; - char *value, *sstr, *dstr, *tstr, *cstr, *nval; - int year, month, day, yday, i, hours, minutes; - //static int mday[12] = {31,28,31,30,31,30,31,31,30,31,30,31}; - int mday[12] = {31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31}; - char value_buffer[VLENGTH + 1]; - - /* Get value and comment from header string */ - value = hgetc(hstring, keyword, value_buffer); - - /* Translate value from ASCII to binary */ - if (value != NULL) - { - sstr = strchr(value, '/'); - dstr = strchr(value, '-'); - - /* Original FITS date format: dd/mm/yy */ - if (sstr > value) - { - *sstr = '\0'; - day = (int) atof(value); - *sstr = '/'; - nval = sstr + 1; - sstr = strchr(nval, '/'); - if (sstr == NULL) - { - sstr = strchr(nval, '-'); - } - if (sstr > value) - { - *sstr = '\0'; - month = (int) atof(nval); - *sstr = '/'; - nval = sstr + 1; - year = (int) atof(nval); - if (day > 31) - { - yday = year; - year = day; - day = yday; - } - if (year >= 0 && year <= 49) - { - year = year + 2000; - } - else if (year < 100) - { - year = year + 1900; - } - if ((year % 4) == 0) - { - mday[1] = 29; - } - else - { - mday[1] = 28; - } - if ((year % 100) == 0 && (year % 400) != 0) - { - mday[1] = 28; - } - if (day > mday[month - 1]) - { - day = mday[month - 1]; - } - else if (day < 1) - { - day = 1; - } - if (mday[1] == 28) - { - yeardays = 365.0; - } - else - { - yeardays = 366.0; - } - yday = day - 1; - for (i = 0; i < month - 1; i++) - { - yday = yday + mday[i]; - } - *dval = (double) year + ((double) yday / yeardays); - return (1); - } - else - { - return (0); - } - } - - /* New FITS date format: yyyy-mm-ddThh:mm:ss[.sss] */ - else if (dstr > value) - { - *dstr = '\0'; - year = (int) atof(value); - *dstr = '-'; - nval = dstr + 1; - dstr = strchr(nval, '-'); - month = 1; - day = 1; - tstr = NULL; - if (dstr > value) - { - *dstr = '\0'; - month = (int) atof(nval); - *dstr = '-'; - nval = dstr + 1; - tstr = strchr(nval, 'T'); - if (tstr > value) - { - *tstr = '\0'; - } - day = (int) atof(nval); - if (tstr > value) - { - *tstr = 'T'; - } - } - - /* If year is < 32, it is really day of month in old format */ - if (year < 32) - { - i = year; - year = day + 1900; - day = i; - } - - if ((year % 4) == 0) - { - mday[1] = 29; - } - else - { - mday[1] = 28; - } - if ((year % 100) == 0 && (year % 400) != 0) - { - mday[1] = 28; - } - if (day > mday[month - 1]) - { - day = mday[month - 1]; - } - else if (day < 1) - { - day = 1; - } - if (mday[1] == 28) - { - yeardays = 365.0; - } - else - { - yeardays = 366.0; - } - yday = day - 1; - for (i = 0; i < month - 1; i++) - { - yday = yday + mday[i]; - } - *dval = (double) year + ((double) yday / yeardays); - - /* Extract time, if it is present */ - if (tstr > value) - { - nval = tstr + 1; - hours = 0.0; - minutes = 0.0; - seconds = 0.0; - cstr = strchr(nval, ':'); - if (cstr > value) - { - *cstr = '\0'; - hours = (int) atof(nval); - *cstr = ':'; - nval = cstr + 1; - cstr = strchr(nval, ':'); - if (cstr > value) - { - *cstr = '\0'; - minutes = (int) atof(nval); - *cstr = ':'; - nval = cstr + 1; - seconds = atof(nval); - } - else - { - minutes = (int) atof(nval); - seconds = 0.0; - } - } - fday = ((3.6e3 * (double) hours) + (6.e1 * (double) minutes) + - seconds) / 8.64e4; - *dval = *dval + (fday / yeardays); - } - return (1); - } - else - { - return (0); - } - } - else - { - return (0); - } -} - - -/* This routine breaks thread safety by using the global var - * "multiline". We likely won't ever call this one, but removed - * it just in case. --PBD - */ -#if 0 -/* Extract IRAF multiple-keyword string value from FITS header string */ - -int -hgetm (hstring, keyword, lstr, str) - -const char *hstring; /* character string containing FITS header information - in the format = {/ } */ -const char *keyword; /* character string containing the root name of the keyword - the value of which is returned. hget searches for a - line beginning with this string. if "[n]" is present, - the n'th token in the value is returned. - (the first 8 characters must be unique) */ -const int lstr; /* Size of str in characters */ -char *str; /* String (returned) */ -{ - char *value; - char *stri; - char keywordi[16]; - int lval, lstri, ikey; - char keyform[8]; - - stri = str; - lstri = lstr; - - sprintf (keywordi, "%s_1", keyword); - if (ksearch (hstring, keywordi)) - strcpy (keyform, "%s_%d"); - else { - sprintf (keywordi, "%s_01", keyword); - if (ksearch (hstring, keywordi)) - strcpy (keyform, "%s_%02d"); - else { - sprintf (keywordi, "%s_001", keyword); - if (ksearch (hstring, keywordi)) - strcpy (keyform, "%s_%03d"); - else - return (0); - } - } - - /* Loop through sequentially-named keywords */ - multiline = 1; - for (ikey = 1; ikey < 500; ikey++) { - char value_buffer[VLENGTH + 1]; - sprintf (keywordi, keyform, keyword, ikey); - - /* Get value for this keyword */ - value = hgetc (hstring, keywordi, value_buffer); - if (value != NULL) { - lval = strlen (value); - if (lval < lstri) - strcpy (stri, value); - else if (lstri > 1) { - strncpy (stri, value, lstri-1); - stri[lstri] = (char) 0; - break; - } - else { - str[0] = value[0]; - break; - } - } - else - break; - stri = stri + lval; - lstri = lstri - lval; - } - multiline = 0; - - /* Return 1 if any keyword found, else 0 */ - if (ikey > 1) - return (1); - else - return (0); -} -#endif - - -/* Extract string value for variable from FITS header string */ - -int -hgetsc(const char *hstring, const char *keyword, const char *wchar, const const int lstr, char *str) - -/* character string containing FITS header information - in the format = {/ } */ -/* character string containing the name of the keyword - the value of which is returned. hget searches for - a line beginning with this string. if "[n]" is - present, the n'th token in the value is returned. - (the first 8 characters must be unique) */ -/* Character of multiple WCS header; =0 if unused */ -/* Size of str in characters */ -/* String (returned) */ -{ - char keyword1[16]; - int lkey; - - if (wchar[0] < (char) 64) - { - return (hgets(hstring, keyword, lstr, str)); - } - else - { - strcpy(keyword1, keyword); - lkey = strlen(keyword); - keyword1[lkey] = wchar[0]; - keyword1[lkey + 1] = (char) 0; - return (hgets(hstring, keyword1, lstr, str)); - } -} -#endif // 0 - -/* Extract string value for variable from FITS header string */ - -static int -hgets(const char *hstring, const char *keyword, const const int lstr, char *str) - -/* character string containing FITS header information -in the format = {/ } */ -/* character string containing the name of the keyword -the value of which is returned. hget searches for a -line beginning with this string. if "[n]" is present, -the n'th token in the value is returned. -(the first 8 characters must be unique) */ -/* Size of str in characters */ -/* String (returned) */ -{ - char *value; - int lval; - char value_buffer[VLENGTH + 1]; - - /* Get value and comment from header string */ - value = hgetc(hstring, keyword, value_buffer); - - if (value != NULL) - { - lval = strlen(value); - if (lval < lstr) - { - strcpy(str, value); - } - else if (lstr > 1) - { - strncpy(str, value, lstr - 1); - } - else - { - str[0] = value[0]; - } - return (1); - } - else - { - return (0); - } -} - -#if 0 -/* Extract number of decimal places for value in FITS header string */ - -int -hgetndec(const char *hstring, const char *keyword, int *ndec) - -/* character string containing FITS header information -in the format = {/ } */ -/* character string containing the name of the keyword -the value of which is returned. hget searches for a -line beginning with this string. if "[n]" is present, -the n'th token in the value is returned. -(the first 8 characters must be unique) */ -/* Number of decimal places in keyword value */ -{ - char *value; - int i, nchar; - char value_buffer[VLENGTH + 1]; - - /* Get value and comment from header string */ - value = hgetc(hstring, keyword, value_buffer); - - /* Find end of string and count backward to decimal point */ - *ndec = 0; - if (value != NULL) - { - nchar = strlen(value); - for (i = nchar - 1; i >= 0; i--) - { - if (value[i] == '.') - { - return (1); - } - *ndec = *ndec + 1; - } - return (1); - } - else - { - return (0); - } -} -#endif // 0 - -/* Extract character value for variable from FITS header string */ - -static char * -hgetc(const char *hstring, const char *keyword0, char *value_buffer) - -/* character string containing FITS header information -in the format = {/ } */ -/* character string containing the name of the keyword -the value of which is returned. hget searches for a -line beginning with this string. if "[n]" is present, -the n'th token in the value is returned. -(the first 8 characters must be unique) */ -{ - //static char cval[80]; - char *cval; - char *value; - char cwhite[2]; - char squot[2], dquot[2], lbracket[2], rbracket[2], slash[2], comma[2]; - char space; - char keyword[81]; /* large for ESO hierarchical keywords */ - char line[100]; - char *vpos, *cpar; - char *q1, *q2, *v1, *v2, *c1, *brack1, *brack2; - //int ipar, i, lkey; - int ipar, i; - -#ifdef USE_SAOLIB - int iel=1, ip=1, nel, np, ier; - char *get_fits_head_str(); - - if( !use_saolib ){ -#endif - if (value_buffer == NULL) - { - return NULL; - } - cval = value_buffer; - - squot[0] = (char) 39; - squot[1] = (char) 0; - dquot[0] = (char) 34; - dquot[1] = (char) 0; - lbracket[0] = (char) 91; - lbracket[1] = (char) 0; - comma[0] = (char) 44; - comma[1] = (char) 0; - rbracket[0] = (char) 93; - rbracket[1] = (char) 0; - slash[0] = (char) 47; - slash[1] = (char) 0; - space = (char) 32; - - /* Find length of variable name */ - strncpy(keyword, keyword0, sizeof(keyword) - 1); - brack1 = strsrch(keyword, lbracket); - if (brack1 == NULL) - { - brack1 = strsrch(keyword, comma); - } - if (brack1 != NULL) - { - *brack1 = '\0'; - brack1++; - } - - /* Search header string for variable name */ - vpos = ksearch(hstring, keyword); - - /* Exit if not found */ - if (vpos == NULL) - { - return (NULL); - } - - /* Initialize line to nulls */ - for (i = 0; i < 100; i++) - { - line[i] = 0; - } - -/* In standard FITS, data lasts until 80th character */ - - /* Extract entry for this variable from the header */ - strncpy(line, vpos, 80); - - /* Check for quoted value */ - q1 = strsrch(line, squot); - c1 = strsrch(line, slash); - if (q1 != NULL) - { - if (c1 != NULL && q1 < c1) - { - q2 = strsrch(q1 + 1, squot); - if (q2 == NULL) - { - q2 = c1 - 1; - while (*q2 == space) - { - q2--; - } - q2++; - } - else if (c1 < q2) - { - c1 = strsrch(q2, slash); - } - } - else if (c1 == NULL) - { - q2 = strsrch(q1 + 1, squot); - if (q2 == NULL) - { - q2 = line + 79; - while (*q2 == space) - { - q2--; - } - q2++; - } - } - else - { - q1 = NULL; - } - } - else - { - q1 = strsrch(line, dquot); - if (q1 != NULL) - { - if (c1 != NULL && q1 < c1) - { - q2 = strsrch(q1 + 1, dquot); - if (q2 == NULL) - { - q2 = c1 - 1; - while (*q2 == space) - { - q2--; - } - q2++; - } - else if (c1 < q2) - { - c1 = strsrch(q2, slash); - } - } - else if (c1 == NULL) - { - q2 = strsrch(q1 + 1, dquot); - if (q2 == NULL) - { - q2 = line + 79; - while (*q2 == space) - { - q2--; - } - q2++; - } - } - else - { - q1 = NULL; - } - } - else - { - q1 = NULL; - q2 = line + 10; - } - } - - /* Extract value and remove excess spaces */ - if (q1 != NULL) - { - v1 = q1 + 1; - v2 = q2; - } - else - { - v1 = strsrch(line, "="); - if (v1 == NULL) - { - v1 = line + 9; - } - else - { - v1 = v1 + 1; - } - c1 = strsrch(line, "/"); - if (c1 != NULL) - { - v2 = c1; - } - else - { - v2 = line + 79; - } - } - - /* Ignore leading spaces if not multiline */ - if (!multiline) - { - while (*v1 == ' ' && v1 < v2) - { - v1++; - } - } - - /* Drop trailing spaces */ - *v2 = '\0'; - if (!multiline) - { - v2--; - while ((*v2 == ' ' || *v2 == (char) 13) && v2 > v1) - { - *v2 = '\0'; - v2--; - } - } - - /* Convert -zero to just plain 0 */ - if (!strcmp(v1, "-0")) - { - v1++; - } - strcpy(cval, v1); - value = cval; - - /* If keyword has brackets, extract appropriate token from value */ - if (brack1 != NULL) - { - brack2 = strsrch(brack1, rbracket); - if (brack2 != NULL) - { - *brack2 = '\0'; - } - if (isnum(brack1)) - { - ipar = atoi(brack1); - cwhite[0] = ' '; - cwhite[1] = '\0'; - if (ipar > 0) - { - for (i = 1; i <= ipar; i++) - { - cpar = strtok(v1, cwhite); - v1 = NULL; - } - if (cpar != NULL) - { - strcpy(cval, cpar); - value = cval; - } - else - { - value = NULL; - } - } - - /* If token counter is negative, include rest of value */ - else if (ipar < 0) - { - for (i = 1; i < -ipar; i++) - { - v1 = strchr(v1, ' '); - if (v1 == NULL) - { - break; - } - else - { - v1 = v1 + 1; - } - } - if (v1 != NULL) - { - strcpy(cval, v1); - value = cval; - } - else - { - value = NULL; - } - } - } - else - { - value = NULL; - /* This seems to be some extended IRAF syntax that - * we won't be using. Removed since it calls igetc. - * --PBD - */ -#if 0 - lkey = strlen (brack1); - for (i = 0; i < lkey; i++) { - if (brack1[i] > 64 && brack1[i] < 91) - brack1[i] = brack1[i] + 32; - } - v1 = igetc (cval, brack1); - if (v1) { - strcpy (cval,v1); - value = cval; - } - else - value = NULL; -#endif - } - } - - return (value); -#ifdef USE_SAOLIB - } else { - return(get_fits_head_str(keyword0, iel, ip, &nel, &np, &ier, hstring)); - } -#endif -} - -#if 0 - -/* Find beginning of fillable blank line before FITS header keyword line */ - -char * -blsearch(const char *hstring, const char *keyword) - -/* Find entry for keyword keyword in FITS header string hstring. - (the keyword may have a maximum of eight letters) - NULL is returned if the keyword is not found */ - -/* character string containing fits-style header -information in the format = {/ } -the default is that each entry is 80 characters long; -however, lines may be of arbitrary length terminated by -nulls, carriage returns or linefeeds, if packed is true. */ -/* character string containing the name of the variable -to be returned. ksearch searches for a line beginning -with this string. The string may be a character -literal or a character variable terminated by a null -or '$'. it is truncated to 8 characters. */ -{ - const char *headlast; - char *loc, *headnext, *pval, *lc, *line; - char *bval; - int icol, nextchar, lkey, nleft, lhstr; - - pval = 0; - - /* Search header string for variable name */ - if (lhead0) - { - lhstr = lhead0; - } - else - { - lhstr = 0; - while (lhstr < 256000 && hstring[lhstr] != 0) - { - lhstr++; - } - } - headlast = hstring + lhstr; - headnext = (char *) hstring; - pval = NULL; - while (headnext < headlast) - { - nleft = headlast - headnext; - loc = strncsrch(headnext, keyword, nleft); - - /* Exit if keyword is not found */ - if (loc == NULL) - { - break; - } - - icol = (loc - hstring) % 80; - lkey = strlen(keyword); - nextchar = (int) *(loc + lkey); - - /* If this is not in the first 8 characters of a line, keep searching */ - if (icol > 7) - { - headnext = loc + 1; - - /* If parameter name in header is longer, keep searching */ - } - else if (nextchar != 61 && nextchar > 32 && nextchar < 127) - { - headnext = loc + 1; - - /* If preceeding characters in line are not blanks, keep searching */ - } - else - { - line = loc - icol; - for (lc = line; lc < loc; lc++) - { - if (*lc != ' ') - { - headnext = loc + 1; - } - } - - /* Return pointer to start of line if match */ - if (loc >= headnext) - { - pval = line; - break; - } - } - } - - /* Return NULL to calling program if keyword is not found */ - if (pval == NULL) - { - return (pval); - } - - /* Return NULL if keyword is found at start of FITS header string */ - if (pval == hstring) - { - return (NULL); - } - - /* Find last nonblank in FITS header string line before requested keyword */ - bval = pval - 80; - while (!strncmp(bval, " ", 8) && bval >= hstring) - { - bval = bval - 80; - } - bval = bval + 80; - - /* Return pointer to calling program if blank lines found */ - if (bval < pval && bval >= hstring) - { - return (bval); - } - else - { - return (NULL); - } -} -#endif // 0 - -/* Find FITS header line containing specified keyword */ - -static char * -ksearch(const char *hstring, const char *keyword) - -/* Find entry for keyword keyword in FITS header string hstring. - (the keyword may have a maximum of eight letters) - NULL is returned if the keyword is not found */ - -/* character string containing fits-style header -information in the format = {/ } -the default is that each entry is 80 characters long; -however, lines may be of arbitrary length terminated by -nulls, carriage returns or linefeeds, if packed is true. */ -/* character string containing the name of the variable -to be returned. ksearch searches for a line beginning -with this string. The string may be a character -literal or a character variable terminated by a null -or '$'. it is truncated to 8 characters. */ -{ - const char *headlast; - char *loc, *headnext, *pval, *lc, *line; - int icol, nextchar, lkey, nleft, lhead, lmax; - -#ifdef USE_SAOLIB - int iel=1, ip=1, nel, np, ier; - char *get_fits_head_str(); - - if( !use_saolib ){ -#endif - - pval = 0; - -/* Find current length of header string */ - if (lhead0) - { - lmax = lhead0; - } - else - { - lmax = 256000; - } - for (lhead = 0; lhead < lmax; lhead++) - { - if (hstring[lhead] == (char) 0) - { - break; - } - } - -/* Search header string for variable name */ - headlast = hstring + lhead; - headnext = (char *) hstring; - pval = NULL; - while (headnext < headlast) - { - nleft = headlast - headnext; - loc = strncsrch(headnext, keyword, nleft); - - /* Exit if keyword is not found */ - if (loc == NULL) - { - break; - } - - icol = (loc - hstring) % 80; - lkey = strlen(keyword); - nextchar = (int) *(loc + lkey); - - /* If this is not in the first 8 characters of a line, keep searching */ - if (icol > 7) - { - headnext = loc + 1; - - /* If parameter name in header is longer, keep searching */ - } - else if (nextchar != 61 && nextchar > 32 && nextchar < 127) - { - headnext = loc + 1; - - /* If preceeding characters in line are not blanks, keep searching */ - } - else - { - line = loc - icol; - for (lc = line; lc < loc; lc++) - { - if (*lc != ' ') - { - headnext = loc + 1; - } - } - - /* Return pointer to start of line if match */ - if (loc >= headnext) - { - pval = line; - break; - } - } - } - -/* Return pointer to calling program */ - return (pval); - -#ifdef USE_SAOLIB - } - else { - if (get_fits_head_str(keyword,iel,ip,&nel,&np,&ier,hstring) != NULL) - return(hstring); - else - return(NULL); - } -#endif -} - -#if 0 - - -/* Return the right ascension in degrees from sexagesimal hours or decimal degrees */ - -double -str2ra(const char *in) - -/* Character string of sexigesimal hours or decimal degrees */ - -{ - double ra; /* Right ascension in degrees (returned) */ - - ra = str2dec(in); - if (strsrch(in, ":")) - { - ra = ra * 15.0; - } - - return (ra); -} - - -/* Return the declination in degrees from sexagesimal or decimal degrees */ - -double -str2dec(const char *in) - -/* Character string of sexigesimal or decimal degrees */ - -{ - double dec; /* Declination in degrees (returned) */ - double deg, min, sec, sign; - char *value, *c1, *c2; - int lval; - char *dchar; - - dec = 0.0; - - /* Return 0.0 if string is null */ - if (in == NULL) - { - return (dec); - } - - /* Translate value from ASCII colon-delimited string to binary */ - if (in[0]) - { - value = (char *) in; - - /* Remove leading spaces */ - while (*value == ' ') - { - value++; - } - - /* Save sign */ - if (*value == '-') - { - sign = -1.0; - value++; - } - else if (*value == '+') - { - sign = 1.0; - value++; - } - else - { - sign = 1.0; - } - - /* Remove trailing spaces */ - lval = strlen(value); - while (value[lval - 1] == ' ') - { - lval--; - } - - if ((c1 = strsrch(value, ":")) == NULL) - { - c1 = strnsrch(value, " ", lval); - } - if (c1 != NULL) - { - *c1 = 0; - deg = (double) atoi(value); - *c1 = ':'; - value = c1 + 1; - if ((c2 = strsrch(value, ":")) == NULL) - { - c2 = strsrch(value, " "); - } - if (c2 != NULL) - { - *c2 = 0; - min = (double) atoi(value); - *c2 = ':'; - value = c2 + 1; - sec = atof(value); - } - else - { - sec = 0.0; - if ((c1 = strsrch(value, ".")) != NULL) - { - min = atof(value); - } - if (strlen(value) > 0) - { - min = (double) atoi(value); - } - } - dec = sign * (deg + (min / 60.0) + (sec / 3600.0)); - } - else if (isnum(value) == 2) - { - if ((dchar = strchr(value, 'D'))) - { - *dchar = 'e'; - } - if ((dchar = strchr(value, 'd'))) - { - *dchar = 'e'; - } - if ((dchar = strchr(value, 'E'))) - { - *dchar = 'e'; - } - dec = sign * atof(value); - } - else - { - dec = sign * (double) atoi(value); - } - } - return (dec); -} -#endif // 0 - - -/* Find string s2 within null-terminated string s1 */ - -static char * -strsrch(const char *s1, const char *s2) - -/* String to search */ -/* String to look for */ - -{ - int ls1; - ls1 = strlen(s1); - return (strnsrch(s1, s2, ls1)); -} - - -/* Find string s2 within string s1 */ - -static char * -strnsrch(const char *s1, const char *s2, const const int ls1) - -/* String to search */ -/* String to look for */ -/* Length of string being searched */ - -{ - char *s, *s1e; - char cfirst, clast; - int i, ls2; - - /* Return null string if either pointer is NULL */ - if (s1 == NULL || s2 == NULL) - { - return (NULL); - } - - /* A zero-length pattern is found in any string */ - ls2 = strlen(s2); - if (ls2 == 0) - { - return ((char *) s1); - } - - /* Only a zero-length string can be found in a zero-length string */ - if (ls1 == 0) - { - return (NULL); - } - - cfirst = (char) s2[0]; - clast = (char) s2[ls2 - 1]; - s1e = (char *) s1 + (int) ls1 - ls2 + 1; - s = (char *) s1; - while (s < s1e) - { - - /* Search for first character in pattern string */ - if (*s == cfirst) - { - - /* If single character search, return */ - if (ls2 == 1) - { - return (s); - } - - /* Search for last character in pattern string if first found */ - if (s[ls2 - 1] == clast) - { - - /* If two-character search, return */ - if (ls2 == 2) - { - return (s); - } - - /* If 3 or more characters, check for rest of search string */ - i = 1; - while (i < ls2 && s[i] == s2[i]) - { - i++; - } - - /* If entire string matches, return */ - if (i >= ls2) - { - return (s); - } - } - } - s++; - } - return (NULL); -} - - -/* Find string s2 within null-terminated string s1 (case-free search) */ - -static char * -strcsrch(const char *s1, const char *s2) - -/* String to search */ -/* String to look for */ - -{ - int ls1; - ls1 = strlen((char *) s1); - return (strncsrch(s1, s2, ls1)); -} - - -/* Find string s2 within string s1 (case-free search) */ - -static char * -strncsrch(const char *s1, const char *s2, const const int ls1) - -/* String to search */ -/* String to look for */ -/* Length of string being searched */ - -{ - char *s, *s1e, sl, *os2; - char cfirst, ocfirst; - char clast = ' '; - char oclast = ' '; - int i, ls2; - - /* Return null string if either pointer is NULL */ - if (s1 == NULL || s2 == NULL) - { - return (NULL); - } - - /* A zero-length pattern is found in any string */ - ls2 = strlen(s2); - if (ls2 == 0) - { - return ((char *) s1); - } - - /* Only a zero-length string can be found in a zero-length string */ - os2 = NULL; - if (ls1 == 0) - { - return (NULL); - } - - /* For one or two characters, set opposite case first and last letters */ - if (ls2 < 3) - { - cfirst = (char) s2[0]; - if (cfirst > 96 && cfirst < 123) - { - ocfirst = cfirst - 32; - } - else if (cfirst > 64 && cfirst < 91) - { - ocfirst = cfirst + 32; - } - else - { - ocfirst = cfirst; - } - if (ls2 > 1) - { - clast = s2[1]; - if (clast > 96 && clast < 123) - { - oclast = clast - 32; - } - else if (clast > 64 && clast < 91) - { - oclast = clast + 32; - } - else - { - oclast = clast; - } - } - } - - /* Else duplicate string with opposite case letters for comparison */ - else - { - os2 = (char *) calloc(ls2, 1); - for (i = 0; i < ls2; i++) - { - if (s2[i] > 96 && s2[i] < 123) - { - os2[i] = s2[i] - 32; - } - else if (s2[i] > 64 && s2[i] < 91) - { - os2[i] = s2[i] + 32; - } - else - { - os2[i] = s2[i]; - } - } - cfirst = s2[0]; - ocfirst = os2[0]; - clast = s2[ls2 - 1]; - oclast = os2[ls2 - 1]; - } - - /* Loop through input string, character by character */ - s = (char *) s1; - s1e = s + (int) ls1 - ls2 + 1; - while (s < s1e) - { - - /* Search for first character in pattern string */ - if (*s == cfirst || *s == ocfirst) - { - - /* If single character search, return */ - if (ls2 == 1) - { - return (s); - } - - /* Search for last character in pattern string if first found */ - sl = s[ls2 - 1]; - if (sl == clast || sl == oclast) - { - - /* If two-character search, return */ - if (ls2 == 2) - { - return (s); - } - - /* If 3 or more characters, check for rest of search string */ - i = 1; - while (i < ls2 && (s[i] == (char) s2[i] || s[i] == os2[i])) - { - i++; - } - - /* If entire string matches, return */ - if (i >= ls2) - { - free(os2); - return (s); - } - } - } - s++; - } - if (os2 != NULL) - { - free(os2); - } - return (NULL); -} - -#if 0 -int -notnum(const char *string) - -/* Character string */ -{ - if (isnum(string)) - { - return (0); - } - else - { - return (1); - } -} -#endif // 0 - - -/* ISNUM-- Return 1 if string is an integer number, - 2 if floating point, - 3 if sexigesimal, with or without decimal point - else 0 - */ - -static int -isnum(const char *string) - -/* Character string */ -{ - int lstr, i, nd, cl; - char cstr, cstr1; - int fpcode; - - /* Return 0 if string is NULL */ - if (string == NULL) - { - return (0); - } - - lstr = strlen(string); - nd = 0; - cl = 0; - fpcode = 1; - - /* Return 0 if string starts with a D or E */ - cstr = string[0]; - if (cstr == 'D' || cstr == 'd' || - cstr == 'E' || cstr == 'e') - { - return (0); - } - - /* Remove trailing spaces */ - while (string[lstr - 1] == ' ') - { - lstr--; - } - - /* Numeric strings contain 0123456789-+ and d or e for exponents */ - for (i = 0; i < lstr; i++) - { - cstr = string[i]; - if (cstr == '\n') - { - break; - } - - /* Ignore leading spaces */ - if (cstr == ' ' && nd == 0) - { - continue; - } - - if ((cstr < 48 || cstr > 57) && - cstr != '+' && cstr != '-' && - cstr != 'D' && cstr != 'd' && - cstr != 'E' && cstr != 'e' && - cstr != ':' && cstr != '.') - { - return (0); - } - else if (cstr == '+' || cstr == '-') - { - if (string[i + 1] == '-' || string[i + 1] == '+') - { - return (0); - } - else if (i > 0) - { - cstr1 = string[i - 1]; - if (cstr1 != 'D' && cstr1 != 'd' && - cstr1 != 'E' && cstr1 != 'e' && - cstr1 != ':' && cstr1 != ' ') - { - return (0); - } - } - } - else if (cstr >= 47 && cstr <= 57) - { - nd++; - - /* Check for colon */ - } - else if (cstr == 58) - { - cl++; - } - if (cstr == '.' || cstr == 'd' || cstr == 'e' || cstr == 'd' || cstr == 'e') - { - fpcode = 2; - } - } - if (nd > 0) - { - if (cl) - { - fpcode = 3; - } - return (fpcode); - } - else - { - return (0); - } -} - -#if 0 - - -/* NUMDEC -- Return number of decimal places in numeric string (-1 if not number) */ - -int -numdec(const char *string) - -/* Numeric string */ -{ - char *cdot; - int lstr; - - if (notnum(string) && !strchr(string, ':')) - { - return (-1); - } - else - { - lstr = strlen(string); - if ((cdot = strchr(string, '.')) == NULL) - { - return (0); - } - else - { - return (lstr - (cdot - string) - 1); - } - } -} - - -#ifdef USE_SAOLIB -int set_saolib(hstring) - void *hstring; -{ - if( *((int *)hstring) == 142857 ) - use_saolib = 1; - else - use_saolib = 0; -} - -#endif - - -/* Remove exponent, leading #, and/or trailing zeroes, if reasonable */ -void -strfix(char *string, int fillblank, int dropzero) - -/* String to modify */ -/* If nonzero, fill blanks with underscores */ -/* If nonzero, drop trailing zeroes */ -{ - char *sdot, *s, *strend, *str, ctemp, *slast; - int ndek, lstr, i; - - /* If number, ignore leading # and remove trailing non-numeric character */ - if (string[0] == '#') - { - strend = string + strlen(string); - str = string + 1; - strend = str + strlen(str) - 1; - ctemp = *strend; - if (!isnum(strend)) - { - *strend = (char) 0; - } - if (isnum(str)) - { - strend = string + strlen(string); - for (str = string; str < strend; str++) - { - *str = *(str + 1); - } - } - else - { - *strend = ctemp; - } - } - - /* Remove positive exponent if there are enough digits given */ - if (isnum(string) > 1 && strsrch(string, "E+") != NULL) - { - lstr = strlen(string); - ndek = (int) (string[lstr - 1] - 48); - ndek = ndek + (10 * ((int) (string[lstr - 2] - 48))); - if (ndek < lstr - 7) - { - lstr = lstr - 4; - string[lstr] = (char) 0; - string[lstr + 1] = (char) 0; - string[lstr + 2] = (char) 0; - string[lstr + 3] = (char) 0; - sdot = strchr(string, '.'); - if (ndek > 0 && sdot != NULL) - { - for (i = 1; i <= ndek; i++) - { - *sdot = *(sdot + 1); - sdot++; - *sdot = '.'; - } - } - } - } - - /* Remove trailing zeroes if they are not significant */ - if (dropzero) - { - if (isnum(string) > 1 && strchr(string, '.') != NULL && - strsrch(string, "E-") == NULL && - strsrch(string, "E+") == NULL && - strsrch(string, "e-") == NULL && - strsrch(string, "e+") == NULL) - { - lstr = strlen(string); - s = string + lstr - 1; - while (*s == '0' && lstr > 1) - { - if (*(s - 1) != '.') - { - *s = (char) 0; - lstr--; - } - s--; - } - } - } - - /* Remove trailing decimal point */ - lstr = strlen(string); - s = string + lstr - 1; - if (*s == '.') - { - *s = (char) 0; - } - - /* Replace embedded blanks with underscores, if requested to */ - if (fillblank) - { - lstr = strlen(string); - slast = string + lstr; - for (s = string; s < slast; s++) - { - if (*s == ' ') - { - *s = '_'; - } - } - } - - return; - -} - -/* Oct 28 1994 New program - * - * Mar 1 1995 Search for / after second quote, not first one - * May 2 1995 Initialize line in HGETC; deal with logicals in HGETL better - * May 4 1995 Declare STRSRCH in KSEARCH - * Aug 7 1995 Fix line initialization in HGETC - * Dec 22 1995 Add HGETRA and HGETDEC to get degrees from xx:xx:xx.xxx string - * - * Jan 26 1996 Fix HGETL to not crash when parameter is not present - * Feb 1 1996 Fix HGETC to deal with quotes correctly - * Feb 1 1996 Fix HGETDEG to deal with sign correctly - * Feb 6 1996 Add HGETS to update character strings - * Feb 8 1996 Fix STRSRCH to find final characters in string - * Feb 23 1996 Add string to degree conversions - * Apr 26 1996 Add HGETDATE to get fractional year from date string - * May 22 1996 Fix documentation; return double from STR2RA and STR2DEC - * May 28 1996 Fix string translation of RA and Dec when no seconds - * Jun 10 1996 Remove unused variables after running lint - * Jun 17 1996 Fix bug which failed to return single character strings - * Jul 1 1996 Skip sign when reading declination after testing for it - * Jul 19 1996 Do not divide by 15 if RA header value is already in degrees - * Aug 5 1996 Add STRNSRCH to search strings which are not null-terminated - * Aug 6 1996 Make minor changes after lint - * Aug 8 1996 Fix ksearch bug which finds wrong keywords - * Aug 13 1996 Fix sign bug in STR2DEC for degrees - * Aug 26 1996 Drop unused variables ICOL0, NLINE, PREVCHAR from KSEARCH - * Sep 10 1996 Fix header length setting code - * Oct 15 1996 Clean up loops and fix ICOL assignment - * Nov 13 1996 Handle integer degrees correctly in STR2DEC - * Nov 21 1996 Make changes for Linux thanks to Sidik Isani - * Dec 12 1996 Add ISNUM to check to see whether strings are numbers - * - * Jan 22 1997 Add ifdefs for Eric Mandel (SAOtng) - * Jan 27 1997 Convert to integer through ATOF so exponents are recognized - * Jul 25 1997 Implement FITS version of ISO date format - * - * Feb 24 1998 Implement code to return IRAF multiple-keyword strings - * Mar 12 1998 Add subroutine NOTNUM - * Mar 27 1998 Add changes to match SKYCAT version - * Apr 30 1998 Add BLSEARCH() to find blank lines before END - * May 27 1998 Add HGETNDEC() to get number of decimal places in entry - * Jun 1 1998 Add VMS patch from Harry Payne at StSci - * Jun 18 1998 Fix code which extracts tokens from string values - * Jul 21 1998 Drop minus sign for values of -0 - * Sep 29 1998 Treat hyphen-separated date as old format if 2-digit year - * Oct 7 1998 Clean up search for last blank line - * - * Apr 5 1999 Check lengths of strings before copying them - * May 5 1999 values.h -> POSIX limits.h: MAXINT->INT_MAX, MAXSHORT->SHRT_MAX - * Jul 15 1999 Add hgetm() options of 1- or 2-digit keyword extensions - * Oct 6 1999 Add gethlength() to return header length - * Oct 14 1999 In ksearch(), search only to null not to end of buffer - * Oct 15 1999 Return 1 from hgetndec() if successful - * Oct 20 1999 Drop unused variable after lint (val in hgetndec) - * Dec 3 1999 Fix isnum() to reject strings starting with a d or e - * Dec 20 1999 Update hgetdate() to get minutes and seconds right - * - * Feb 10 2000 Parse RA and Dec with spaces as well as colons as separators - * Feb 11 2000 Add null at end of multi-line keyword value character string - * Feb 25 2000 Change max search string length from 57600 to 256000 - * Mar 15 2000 Deal with missing second quotes in string values - * Mar 17 2000 Return 2 from isnum() if number is floating point (.de) - * Mar 17 2000 Ignore leading # for numeric values in header - * Mar 21 2000 Implement -n to get string value starting with nth token - * Apr 5 2000 Reject +- in isnum() - * Jun 9 2000 Read keyword values even if no equal sign is present - * Sep 20 2000 Ignore linefeed at end of number in isnum() - * Oct 23 2000 Fix handling of embedded + or - in isnum() - * - * Jan 19 2000 Return 0 from isnum(), str2ra(), and str2dec() if string is null - * Mar 30 2001 Fix header length finding algorithm in ksearch() - * Jul 13 2001 Make val[] static int instead of int; drop unused variables - * Sep 12 2001 Read yyyy/mm/dd dates as well as dd/mm/yyyy - * Sep 20 2001 Ignore leading spaces in str2dec() - * Sep 20 2001 Ignore trailing spaces in isnum() - * - * Apr 3 2002 Add hgetr8c(), hgeti4c(), and hgetsc() for multiple WCS handling - * Apr 26 2002 Fix bug in hgetsc(), hgeti4c(), and hgetr8c() found by Bill Joye - * Jun 26 2002 Do not drop leading or trailing spaces in multi-line values - * Aug 6 2002 Add strcsrch() and strncsrch() for case-insensitive searches - * Aug 30 2002 Fix bug so strcsrch() really is case-insensitive - * Oct 20 2003 Add numdec() to return number of decimal places in a string - * Dec 9 2003 Fix numdec() to return 0 if no digits after decimal point - * - * Feb 26 2004 Extract value from keyword=value strings within a keyword value - * Apr 9 2004 Use strncsrch() in ksearch() to find differently-cased keywords - * Apr 28 2004 Free os2 in strncsrch() only if it is allocated - * Jul 13 2004 Accept D, d, E, or e as exponent delimiter in floating points - * Aug 30 2004 Change numdec() to accept sexigesimal numbers (:'s) - * - * Jun 27 2005 Drop unused variables - * Aug 30 2005 Adjust code in hlength() - * - * Jun 20 2006 Initialize uninitialized variables in strnsrch() - * Jun 29 2006 Add new subroutine strfix() to clean strings for other uses - * Jul 13 2006 Increase maximum number of multiline keywords from 20 to 500 - * - * Jan 4 2007 Declare header, keyword to be const - * Jan 4 2007 Change WCS letter from char to char* - * Feb 28 2007 If header length is not set in hlength, set it to 0 - * May 31 2007 Add return value of 3 to isnum() if string has colon(s) - * Aug 22 2007 If closing quote not found, make one up - * Sep 6 2016 Added third arg to hgetc() to correct a 'return ptr to stack' issue. - */ -#endif // 0 - -#endif // HGET_H diff --git a/rawspec_rawutils.c b/rawspec_rawutils.c index f22e948..b0a0d64 100644 --- a/rawspec_rawutils.c +++ b/rawspec_rawutils.c @@ -5,7 +5,6 @@ #include #include "rawspec_rawutils.h" -#include "hget.h" int32_t rawspec_raw_get_s32(const char * buf, const char * key, int32_t def) { @@ -151,38 +150,89 @@ int rawspec_raw_header_size(char * hdr, size_t len, int directio) return 0; } +const uint64_t KEY_UINT64_BLOCSIZE = GUPPI_RAW_KEY_UINT64_ID_LE('B','L','O','C','S','I','Z','E'); +const uint64_t KEY_UINT64_NPOL = GUPPI_RAW_KEY_UINT64_ID_LE('N','P','O','L',' ',' ',' ',' '); +const uint64_t KEY_UINT64_OBSNCHAN = GUPPI_RAW_KEY_UINT64_ID_LE('O','B','S','N','C','H','A','N'); +const uint64_t KEY_UINT64_NBITS = GUPPI_RAW_KEY_UINT64_ID_LE('N','B','I','T','S',' ',' ',' '); +const uint64_t KEY_UINT64_OBSFREQ = GUPPI_RAW_KEY_UINT64_ID_LE('O','B','S','F','R','E','Q',' '); +const uint64_t KEY_UINT64_OBSBW = GUPPI_RAW_KEY_UINT64_ID_LE('O','B','S','B','W',' ',' ',' '); +const uint64_t KEY_UINT64_TBIN = GUPPI_RAW_KEY_UINT64_ID_LE('T','B','I','N',' ',' ',' ',' '); +const uint64_t KEY_UINT64_DIRECTIO = GUPPI_RAW_KEY_UINT64_ID_LE('D','I','R','E','C','T','I','O'); +const uint64_t KEY_UINT64_PKTIDX = GUPPI_RAW_KEY_UINT64_ID_LE('P','K','T','I','D','X',' ',' '); +const uint64_t KEY_UINT64_BEAM_ID = GUPPI_RAW_KEY_UINT64_ID_LE('B','E','A','M','_','I','D',' '); +const uint64_t KEY_UINT64_NBEAM = GUPPI_RAW_KEY_UINT64_ID_LE('N','B','E','A','M',' ',' ',' '); +const uint64_t KEY_UINT64_NANTS = GUPPI_RAW_KEY_UINT64_ID_LE('N','A','N','T','S',' ',' ',' '); +const uint64_t KEY_UINT64_RA_STR = GUPPI_RAW_KEY_UINT64_ID_LE('R','A','_','S','T','R',' ',' '); +const uint64_t KEY_UINT64_DEC_STR = GUPPI_RAW_KEY_UINT64_ID_LE('D','E','C','_','S','T','R',' '); +const uint64_t KEY_UINT64_STT_IMJD = GUPPI_RAW_KEY_UINT64_ID_LE('S','T','T','_','I','M','J','D'); +const uint64_t KEY_UINT64_STT_SMJD = GUPPI_RAW_KEY_UINT64_ID_LE('S','T','T','_','S','M','J','D'); +const uint64_t KEY_UINT64_SRC_NAME = GUPPI_RAW_KEY_UINT64_ID_LE('S','R','C','_','N','A','M','E'); +const uint64_t KEY_UINT64_TELESCOP = GUPPI_RAW_KEY_UINT64_ID_LE('T','E','L','E','S','C','O','P'); + +void _rawspec_header_parse_metadata(const char* entry, void* raw_hdr_void) { + rawspec_raw_hdr_t * raw_hdr = (rawspec_raw_hdr_t*) raw_hdr_void; + + if(((uint64_t*)entry)[0] == KEY_UINT64_BLOCSIZE) + hgeti4(entry, "BLOCSIZE", &raw_hdr->blocsize); + else if(((uint64_t*)entry)[0] == KEY_UINT64_NPOL) + hgeti4(entry, "NPOL", &raw_hdr->npol); + else if(((uint64_t*)entry)[0] == KEY_UINT64_OBSNCHAN) + hgeti4(entry, "OBSNCHAN", &raw_hdr->obsnchan); + else if(((uint64_t*)entry)[0] == KEY_UINT64_NBITS) + hgetu4(entry, "NBITS", &raw_hdr->nbits); + else if(((uint64_t*)entry)[0] == KEY_UINT64_OBSFREQ) + hgetr8(entry, "OBSFREQ", &raw_hdr->obsfreq); + else if(((uint64_t*)entry)[0] == KEY_UINT64_OBSBW) + hgetr8(entry, "OBSBW", &raw_hdr->obsbw); + else if(((uint64_t*)entry)[0] == KEY_UINT64_TBIN) + hgetr8(entry, "TBIN", &raw_hdr->tbin); + else if(((uint64_t*)entry)[0] == KEY_UINT64_DIRECTIO) + hgeti4(entry, "DIRECTIO", &raw_hdr->directio); + else if(((uint64_t*)entry)[0] == KEY_UINT64_PKTIDX) + hgetu8(entry, "PKTIDX", &raw_hdr->pktidx); + else if(((uint64_t*)entry)[0] == KEY_UINT64_BEAM_ID) + hgeti4(entry, "BEAM_ID", &raw_hdr->beam_id); + else if(((uint64_t*)entry)[0] == KEY_UINT64_NBEAM) + hgeti4(entry, "NBEAM", &raw_hdr->nbeam); + else if(((uint64_t*)entry)[0] == KEY_UINT64_NANTS) + hgetu4(entry, "NANTS", &raw_hdr->nants); + else if(((uint64_t*)entry)[0] == KEY_UINT64_RA_STR) { + char tmp[72]; + hgets(entry, "RA_STR", 72, tmp); + raw_hdr->ra = rawspec_raw_hmsstr_to_h(tmp); + } + else if(((uint64_t*)entry)[0] == KEY_UINT64_DEC_STR) { + char tmp[72]; + hgets(entry, "DEC_STR", 72, tmp); + raw_hdr->dec = rawspec_raw_dmsstr_to_d(tmp); + } + else if(((uint64_t*)entry)[0] == KEY_UINT64_STT_IMJD) { + double tmp = 0.0; + hgetr8(entry, "STT_IMJD", &tmp); + raw_hdr->mjd += tmp/86400.0; + } + else if(((uint64_t*)entry)[0] == KEY_UINT64_STT_SMJD) { + double tmp = 0.0; + hgetr8(entry, "STT_SMJD", &tmp); + raw_hdr->mjd += tmp/86400.0; + } + else if(((uint64_t*)entry)[0] == KEY_UINT64_SRC_NAME) + hgets(entry, "SRC_NAME", 72, raw_hdr->src_name); + else if(((uint64_t*)entry)[0] == KEY_UINT64_TELESCOP) + hgets(entry, "TELESCOP", 72, raw_hdr->telescop); +} + // Parses rawspec related RAW header params from buf into raw_hdr. void rawspec_raw_parse_header(const char * buf, rawspec_raw_hdr_t * raw_hdr) { - int smjd; - int imjd; - char tmp[80]; - - raw_hdr->blocsize = rawspec_raw_get_s32(buf, "BLOCSIZE", 0); - raw_hdr->npol = rawspec_raw_get_s32(buf, "NPOL", 0); - raw_hdr->obsnchan = rawspec_raw_get_s32(buf, "OBSNCHAN", 0); - raw_hdr->nbits = rawspec_raw_get_u32(buf, "NBITS", 8); - raw_hdr->obsfreq = rawspec_raw_get_dbl(buf, "OBSFREQ", 0.0); - raw_hdr->obsbw = rawspec_raw_get_dbl(buf, "OBSBW", 0.0); - raw_hdr->tbin = rawspec_raw_get_dbl(buf, "TBIN", 0.0); - raw_hdr->directio = rawspec_raw_get_s32(buf, "DIRECTIO", 0); - raw_hdr->pktidx = rawspec_raw_get_u64(buf, "PKTIDX", -1); - raw_hdr->beam_id = rawspec_raw_get_s32(buf, "BEAM_ID", -1); - raw_hdr->nbeam = rawspec_raw_get_s32(buf, "NBEAM", -1); - raw_hdr->nants = rawspec_raw_get_u32(buf, "NANTS", 1); - - rawspec_raw_get_str(buf, "RA_STR", "0.0", tmp, 80); - raw_hdr->ra = rawspec_raw_hmsstr_to_h(tmp); - - rawspec_raw_get_str(buf, "DEC_STR", "0.0", tmp, 80); - raw_hdr->dec = rawspec_raw_dmsstr_to_d(tmp); - - imjd = rawspec_raw_get_s32(buf, "STT_IMJD", 51545); // TODO use double? - smjd = rawspec_raw_get_s32(buf, "STT_SMJD", 0); // TODO use double? - raw_hdr->mjd = ((double)imjd) + ((double)smjd)/86400.0; - - rawspec_raw_get_str(buf, "SRC_NAME", "Unknown", raw_hdr->src_name, 80); - rawspec_raw_get_str(buf, "TELESCOP", "Unknown", raw_hdr->telescop, 80); + snprintf(raw_hdr->src_name, 72, "Unknown"); + snprintf(raw_hdr->telescop, 72, "Unknown"); + raw_hdr->mjd = 0.0; + + guppiraw_metadata_t metadata = {0}; + metadata.user_data = raw_hdr; + metadata.user_callback = _rawspec_header_parse_metadata; + guppiraw_header_string_parse_metadata(&metadata, buf, -1); } // Reads RAW file params from fd. On entry, fd is assumed to be at the start diff --git a/rawspec_rawutils.h b/rawspec_rawutils.h index 5bf1a1e..b48b960 100644 --- a/rawspec_rawutils.h +++ b/rawspec_rawutils.h @@ -4,6 +4,8 @@ #include #include +#include "guppirawc99/header.h" + typedef struct { int directio; size_t blocsize; From cd8a48d038b77fb342cc54bf84fd7cf2295e9e74 Mon Sep 17 00:00:00 2001 From: radonnachie Date: Fri, 29 Jul 2022 16:05:31 +0000 Subject: [PATCH 3/9] +* rawspec_raw_read_guppiraw_header --- rawspec.c | 6 ++++-- rawspec_rawutils.c | 27 ++++++++++++++++++++++++--- rawspec_rawutils.h | 1 + 3 files changed, 29 insertions(+), 5 deletions(-) diff --git a/rawspec.c b/rawspec.c index 612a214..77cb057 100644 --- a/rawspec.c +++ b/rawspec.c @@ -200,6 +200,8 @@ int main(int argc, char *argv[]) size_t total_bytes_read; off_t pos; rawspec_raw_hdr_t raw_hdr; + guppiraw_header_t guppiraw_header; + guppiraw_header.metadata.user_data = &raw_hdr; callback_data_t cb_data[MAX_OUTPUTS]; rawspec_context ctx; int ant = -1; @@ -572,7 +574,7 @@ int main(int argc, char *argv[]) posix_fadvise(fdin, 0, 0, POSIX_FADV_SEQUENTIAL); // Read obs params - pos = rawspec_raw_read_header(fdin, &raw_hdr); + pos = rawspec_raw_read_guppiraw_header(fdin, &guppiraw_header); if(pos <= 0) { if(pos == -1) { fprintf(stderr, "error getting obs params from %s\n", fname); @@ -1031,7 +1033,7 @@ int main(int argc, char *argv[]) bi++; // Read obs params of next block - pos = rawspec_raw_read_header(fdin, &raw_hdr); + pos = rawspec_raw_read_guppiraw_header(fdin, &guppiraw_header); if(pos <= 0) { if(pos == -1) { fprintf(stderr, "error getting obs params from %s [%s]\n", diff --git a/rawspec_rawutils.c b/rawspec_rawutils.c index b0a0d64..3b3f332 100644 --- a/rawspec_rawutils.c +++ b/rawspec_rawutils.c @@ -240,8 +240,7 @@ void rawspec_raw_parse_header(const char * buf, rawspec_raw_hdr_t * raw_hdr) // of the subsequent data block and the file descriptor `fd` will also refer to // that location in the file. On EOF, this function returns 0. On failure, // this function returns -1 and the location to which fd refers is undefined. -off_t rawspec_raw_read_header(int fd, rawspec_raw_hdr_t * raw_hdr) -{ +off_t rawspec_raw_read_guppiraw_header(int fd, guppiraw_header_t * header) { int i; // Ensure that hdr is aligned to a 512-byte boundary so that it can be used // with files opened with O_DIRECT. @@ -258,7 +257,17 @@ off_t rawspec_raw_read_header(int fd, rawspec_raw_hdr_t * raw_hdr) return 0; } - rawspec_raw_parse_header(hdr, raw_hdr); + if(header->metadata.user_data == NULL) { + header->metadata.user_data = malloc(sizeof(rawspec_raw_hdr_t)); + } + header->metadata.user_callback = _rawspec_header_parse_metadata; + + rawspec_raw_hdr_t* raw_hdr = (rawspec_raw_hdr_t*) header->metadata.user_data; + snprintf(raw_hdr->src_name, 72, "Unknown"); + snprintf(raw_hdr->telescop, 72, "Unknown"); + raw_hdr->mjd = 0.0; + + guppiraw_header_parse(header, hdr, hdr_size); if(raw_hdr->blocsize == 0) { fprintf(stderr, " BLOCSIZE not found in header\n"); @@ -308,3 +317,15 @@ off_t rawspec_raw_read_header(int fd, rawspec_raw_hdr_t * raw_hdr) return pos; } + +// Reads RAW file params from fd. On entry, fd is assumed to be at the start +// of a RAW header section. On success, this function returns the file offset +// of the subsequent data block and the file descriptor `fd` will also refer to +// that location in the file. On EOF, this function returns 0. On failure, +// this function returns -1 and the location to which fd refers is undefined. +off_t rawspec_raw_read_header(int fd, rawspec_raw_hdr_t * raw_hdr) +{ + guppiraw_header_t header = {0}; + header.metadata.user_data = (void*) raw_hdr; + return rawspec_raw_read_guppiraw_header(fd, &header); +} diff --git a/rawspec_rawutils.h b/rawspec_rawutils.h index b48b960..0e73673 100644 --- a/rawspec_rawutils.h +++ b/rawspec_rawutils.h @@ -63,6 +63,7 @@ void rawspec_raw_parse_header(const char * buf, rawspec_raw_hdr_t * raw_hdr); // that location in the file. On EOF, this function returns 0. On failure, // this function returns -1 and the location to which fd refers is undefined. off_t rawspec_raw_read_header(int fd, rawspec_raw_hdr_t * raw_hdr); +off_t rawspec_raw_read_guppiraw_header(int fd, guppiraw_header_t * header); #ifdef __cplusplus } From 7466c64b69043764b79df1073905e5921a0b3e88 Mon Sep 17 00:00:00 2001 From: radonnachie Date: Sat, 30 Jul 2022 09:27:50 +0000 Subject: [PATCH 4/9] ++ write FILE_FORMAT_GUPPIRAW --- rawspec.c | 49 ++++++++++++++++++++++++-- rawspec_callback.h | 3 ++ rawspec_file.c | 85 ++++++++++++++++++++++++++++++++++++++-------- 3 files changed, 120 insertions(+), 17 deletions(-) diff --git a/rawspec.c b/rawspec.c index 77cb057..14254ec 100644 --- a/rawspec.c +++ b/rawspec.c @@ -115,6 +115,7 @@ void usage(const char *argv0) { " -i, --ics=W1[,W2...] Output incoherent-sum (exclusively, unless with -S)\n" " specifying per antenna-weights or a singular, uniform weight\n" " -j, --fbh5 Format output Filterbank files as FBH5 (.h5) instead of SIGPROC(.fil)\n" + " -k, --guppi Output GUPPI RAW files (.raw) instead of Filterbank files\n" " -n, --nchan=N Number of coarse channels to process [all]\n" " -o, --outidx=N First index number for output files [0]\n" " -p --pols={1|4}[,...] Number of output polarizations [1]\n" @@ -239,7 +240,7 @@ int main(int argc, char *argv[]) // Parse command line. argv0 = argv[0]; - while((opt=getopt_long(argc, argv, "a:b:d:f:g:HSjzs:i:n:o:p:r:t:hv", long_opts, NULL)) != -1) { + while((opt=getopt_long(argc, argv, "a:b:d:f:g:HSjkzs:i:n:o:p:r:t:hv", long_opts, NULL)) != -1) { switch (opt) { case 'h': // Help usage(argv0); @@ -250,6 +251,10 @@ int main(int argc, char *argv[]) flag_file_output = FILE_FORMAT_FBH5; break; + case 'k': // GUPPI RAW output format requested + flag_file_output = FILE_FORMAT_GUPPIRAW; + break; + case 'z': // Selected dynamic debugging flag_debugging = 1; break; @@ -412,6 +417,9 @@ int main(int argc, char *argv[]) case FILE_FORMAT_FBSIGPROC: printf("writing output files in SIGPROC Filterbank format\n"); break; + case FILE_FORMAT_GUPPIRAW: + printf("writing GUPPI RAW output files\n"); + break; } } @@ -627,7 +635,7 @@ int main(int argc, char *argv[]) fprintf(stderr, "OBSBW = %g\n", raw_hdr.obsbw); fprintf(stderr, "TBIN = %g\n", raw_hdr.tbin); #endif // VERBOSE - if(raw_hdr.nants > 1 && !(per_ant_out || ctx.incoherently_sum)){ + if(raw_hdr.nants > 1 && !(per_ant_out || ctx.incoherently_sum) && flag_file_output != FILE_FORMAT_GUPPIRAW){ printf("NANTS = %d >1: Enabling --split-ant in lieu of neither --split-ant nor --ics flags.\n", raw_hdr.nants); per_ant_out = 1; } @@ -795,6 +803,43 @@ int main(int argc, char *argv[]) printf("output %d Nds = %u, Nf = %u\n", i, cb_data[i].Nds, cb_data[i].Nf); } cb_data[i].Nant = raw_hdr.nants; + + if(flag_file_output == FILE_FORMAT_GUPPIRAW) { + if(ctx.incoherently_sum) { + guppiraw_header_copy(&cb_data[i].guppiraw_header_ics, &guppiraw_header); + cb_data[i].guppiraw_header_ics.metadata.user_data = malloc(sizeof(rawspec_raw_hdr_t)); + memcpy(cb_data[i].guppiraw_header_ics.metadata.user_data, &raw_hdr, sizeof(rawspec_raw_hdr_t)); + guppiraw_header_put_string(&cb_data[i].guppiraw_header_ics, "DATATYPE", "FLOAT"); + + cb_data[i].guppiraw_header_ics.metadata.datashape.block_size = ctx.h_pwrbuf_size[i]/ctx.Nant; + cb_data[i].guppiraw_header_ics.metadata.datashape.n_ant = 1; + cb_data[i].guppiraw_header.metadata.datashape.n_bit = (sizeof(float)*8)/2; // to account for assumption of data being complex + cb_data[i].guppiraw_header_ics.metadata.datashape.n_time = ctx.Nds[i]; + cb_data[i].guppiraw_header_ics.metadata.datashape.n_pol = ctx.Npolout[i]; + cb_data[i].guppiraw_header_ics.metadata.datashape.n_aspectchan = ctx.Nts[i] * ctx.Nc/ctx.Nant; + guppiraw_header_put_metadata(&cb_data[i].guppiraw_header); + } + + guppiraw_header_copy(&cb_data[i].guppiraw_header, &guppiraw_header); + cb_data[i].guppiraw_header.metadata.user_data = malloc(sizeof(rawspec_raw_hdr_t)); + memcpy(cb_data[i].guppiraw_header.metadata.user_data, &raw_hdr, sizeof(rawspec_raw_hdr_t)); + + if(cb_data[i].per_ant_out) { + cb_data[i].guppiraw_header.metadata.datashape.n_ant = 1; + } + cb_data[i].guppiraw_header.metadata.datashape.n_bit = (sizeof(float)*8)/2; // to account for assumption of data being complex + cb_data[i].guppiraw_header.metadata.datashape.n_time = ctx.Nds[i]; + cb_data[i].guppiraw_header.metadata.datashape.n_pol = ctx.Npolout[i]; + cb_data[i].guppiraw_header.metadata.datashape.n_aspectchan = ctx.Nts[i] * ctx.Nc/ctx.Nant; + guppiraw_header_put_string(&cb_data[i].guppiraw_header, "DATATYPE", "FLOAT"); + guppiraw_header_put_metadata(&cb_data[i].guppiraw_header); + char* header_string = guppiraw_header_malloc_string(&guppiraw_header); + printf("```%s```\n", header_string); + free(header_string); + header_string = guppiraw_header_malloc_string(&cb_data[i].guppiraw_header); + printf("```%s```\n", header_string); + free(header_string); + } } #if 0 if(output_mode == RAWSPEC_NET) { diff --git a/rawspec_callback.h b/rawspec_callback.h index 4295693..52f5d53 100644 --- a/rawspec_callback.h +++ b/rawspec_callback.h @@ -4,6 +4,7 @@ #include #include "hdf5.h" #include "rawspec_fbutils.h" +#include "guppirawc99/header.h" typedef struct { int active; // Still active? 1=yes, 0=no @@ -61,6 +62,8 @@ typedef struct { unsigned int exit_soon; // Added for GUPPI RAW output 2022-07 + guppiraw_header_t guppiraw_header; + guppiraw_header_t guppiraw_header_ics; } callback_data_t; #endif // _RAWSPEC_CALLBACK_H_ diff --git a/rawspec_file.c b/rawspec_file.c index 6966a3e..f9805fa 100644 --- a/rawspec_file.c +++ b/rawspec_file.c @@ -7,6 +7,7 @@ #include "rawspec_file.h" #include "fbh5_defs.h" +#include "guppirawc99.h" // Open a single Filterbank file for one of the following: // * nants = 0 @@ -28,6 +29,9 @@ int open_output_file(callback_data_t *cb_data, const char * dest, const char *st case FILE_FORMAT_FBSIGPROC: strcpy(fileext, "fil"); break; + case FILE_FORMAT_GUPPIRAW: + strcpy(fileext, "raw"); + break; } if(dest && dest[0]) { // Look for last '/' in stem @@ -73,6 +77,15 @@ int open_output_file(callback_data_t *cb_data, const char * dest, const char *st } return fd; + case FILE_FORMAT_GUPPIRAW: + // Open a SIGPROC Filterbank output file. + fd = open(fname, O_WRONLY | O_CREAT | O_TRUNC, 0664); + if(fd == -1) { + cb_data->exit_soon = 1; + perror(fname); + } + return fd; + } } @@ -114,9 +127,41 @@ void * dump_file_thread_func(void *arg) { callback_data_t * cb_data = (callback_data_t *)arg; int retcode; + const size_t spectra_stride = cb_data->h_pwrbuf_size / (cb_data->Nds * sizeof(float)); + const size_t pol_stride = spectra_stride / cb_data->fb_hdr.nifs; + const size_t ant_stride = pol_stride / cb_data->Nant; + const size_t chan_stride = ant_stride / cb_data->Nf; // Single antenna case - if(cb_data->fd && cb_data->h_pwrbuf && (cb_data->Nant == 1)) { + if(cb_data->fd && cb_data->h_pwrbuf) { + if(cb_data->flag_file_output == FILE_FORMAT_GUPPIRAW) { + // printf("guppiraw_write antenna: %ld bytes\n\t", cb_data->h_pwrbuf_size / (cb_data->per_ant_out ? cb_data->Nant : 1 )); + for(size_t i = 0; i < (cb_data->per_ant_out ? cb_data->Nant : 1); i++){ + // printf("%d ", cb_data->fd[i]); + if(cb_data->fd[i] == -1){ + // Assume that the following file-descriptors aren't valid + break; + } + size_t bytes_written = guppiraw_write_block_arbitrary( + cb_data->fd[i], + &cb_data->guppiraw_header, + cb_data->h_pwrbuf + ant_stride*i, + ant_stride*sizeof(float), + chan_stride*sizeof(float), + spectra_stride*sizeof(float), + pol_stride*sizeof(float) + ); + if(bytes_written + < cb_data->h_pwrbuf_size / (cb_data->per_ant_out ? cb_data->Nant : 1) + ) { + cb_data->exit_soon = 1; + cb_data->output_thread_valid = 0; + fprintf(stderr, "GUPPIRAW-WRITE-ERROR: wrote %ld (%ld vs %ld)\n", bytes_written, cb_data->h_pwrbuf_size, cb_data->guppiraw_header.metadata.datashape.block_size); + } + } + // printf("\n"); + } + else if(cb_data->Nant == 1) { if(cb_data->debug_callback) printf("dump_file_thread_func: write for nants=0\n"); @@ -142,15 +187,8 @@ void * dump_file_thread_func(void *arg) } break; } - } - - // Multiple antennas, split output - if(cb_data->fd && cb_data->h_pwrbuf && (cb_data->Nant > 1)) { - if(cb_data->per_ant_out) { - size_t spectra_stride = cb_data->h_pwrbuf_size / (cb_data->Nds * sizeof(float)); - size_t pol_stride = spectra_stride / cb_data->fb_hdr.nifs; - size_t ant_stride = pol_stride / cb_data->Nant; - + } else if(cb_data->per_ant_out && cb_data->Nant > 1) { + // Multiple antennas, split output for(size_t k = 0; k < cb_data->Nds; k++){// Spectra out for(size_t j = 0; j < cb_data->fb_hdr.nifs; j++){// Npolout for(size_t i = 0; i < cb_data->Nant; i++){ @@ -164,10 +202,10 @@ void * dump_file_thread_func(void *arg) switch(cb_data->flag_file_output) { case FILE_FORMAT_FBH5: retcode = fbh5_write(&(cb_data->fbh5_ctx_ant[i]), - &(cb_data->fb_hdr), - cb_data->h_pwrbuf + i * ant_stride + j * pol_stride + k * spectra_stride, - ant_stride * sizeof(float), - cb_data->debug_callback); + &(cb_data->fb_hdr), + cb_data->h_pwrbuf + i * ant_stride + j * pol_stride + k * spectra_stride, + ant_stride * sizeof(float), + cb_data->debug_callback); if(retcode != 0) { cb_data->exit_soon = 1; cb_data->output_thread_valid = 0; @@ -186,7 +224,7 @@ void * dump_file_thread_func(void *arg) } // for(size_t i = 0; i < cb_data->Nant; i++) } // for(size_t j = 0; j < cb_data->fb_hdr.nifs; j++) } // for(size_t k = 0; k < cb_data->Nds; k++) - } // if(cb_data->per_ant_out) + } // else if(cb_data->per_ant_out) } // if(cb_data->fd && cb_data->h_pwrbuf) // Multiple antennas, ICS output @@ -215,6 +253,23 @@ void * dump_file_thread_func(void *arg) fprintf(stderr, "SIGPROC-WRITE-ERROR\n"); } break; + case FILE_FORMAT_GUPPIRAW: + if( + guppiraw_write_block_arbitrary( + cb_data->fd_ics, + &cb_data->guppiraw_header_ics, + cb_data->h_pwrbuf, + ant_stride*sizeof(float), + chan_stride*sizeof(float), + spectra_stride*sizeof(float), + pol_stride*sizeof(float) + ) < cb_data->h_pwrbuf_size / (cb_data->per_ant_out ? 1 : cb_data->Nant) + ) { + cb_data->exit_soon = 1; + cb_data->output_thread_valid = 0; + fprintf(stderr, "GUPPIRAW-WRITE-ERROR\n"); + } + break; } } From 78e7020a9dd6d2bf4034c3f5460e8c223c9c4980 Mon Sep 17 00:00:00 2001 From: radonnachie Date: Mon, 1 Aug 2022 15:04:34 +0000 Subject: [PATCH 5/9] ++ guppiraw complex output, WIP --- rawspec.c | 39 +++++++++---- rawspec.h | 3 + rawspec_file.c | 57 +++++++++++-------- rawspec_gpu.cu | 149 +++++++++++++++++++++++++++++++++++++------------ 4 files changed, 177 insertions(+), 71 deletions(-) diff --git a/rawspec.c b/rawspec.c index 14254ec..4b42daf 100644 --- a/rawspec.c +++ b/rawspec.c @@ -233,7 +233,7 @@ int main(int argc, char *argv[]) // Init rawspec context memset(&ctx, 0, sizeof(ctx)); - ctx.Npolout[0] = 1; // others will be set later + ctx.Npolout[0] = 0; // others will be set later // Exit status after mallocs have occured. int exit_status = 0; @@ -465,21 +465,32 @@ int main(int argc, char *argv[]) ctx.Nas[2] = 3072; } + if(ctx.Npolout[0] == 0) { + // default Npolout (-1 means match input, only for GUPPIRAW) + ctx.Npolout[0] = (flag_file_output == FILE_FORMAT_GUPPIRAW ? -1 : 1); + } + // Validate polout values for(i=0; i 0) { // Copy value from previous output product ctx.Npolout[i] = ctx.Npolout[i-1]; - } else if(ctx.Npolout[i]!=1 && abs(ctx.Npolout[i])!=4) { - fprintf(stderr, - "error: number of output pols must be 1 or +/- 4\n"); - return 1; + } else { + if(flag_file_output == FILE_FORMAT_GUPPIRAW && ctx.Npolout[i]!=-1) { + fprintf(stderr, + "error: GUPPI RAW output mode necessitates auto-output pols (-1)\n"); + return 1; + } else if(abs(ctx.Npolout[i])!=1 && abs(ctx.Npolout[i])!=4) { + fprintf(stderr, + "error: number of output pols must be 1 or +/- 4\n"); + return 1; + } } // Full-pol mode is not supported for network output if(ctx.Npolout[i] != 1 && output_mode != RAWSPEC_FILE) { fprintf(stderr, - "error: full-pol mode is not supported for network output\n"); + "error: full mode is not supported for network output\n"); return 1; } } @@ -639,6 +650,15 @@ int main(int argc, char *argv[]) printf("NANTS = %d >1: Enabling --split-ant in lieu of neither --split-ant nor --ics flags.\n", raw_hdr.nants); per_ant_out = 1; } + if(flag_file_output == FILE_FORMAT_GUPPIRAW) { + // assume ctx.Npolout == -1 due to prior checks + for(i=0; i 2, {2*N + 1 -> 1} + // negative to indicate complex output (ctx.complex_output) + ctx.complex_output = 1; + ctx.Npolout[i] = - (2 - (raw_hdr.npol % 2)); + } + } // If splitting output per antenna, re-alloc the fd array. if(per_ant_out) { @@ -811,11 +831,9 @@ int main(int argc, char *argv[]) memcpy(cb_data[i].guppiraw_header_ics.metadata.user_data, &raw_hdr, sizeof(rawspec_raw_hdr_t)); guppiraw_header_put_string(&cb_data[i].guppiraw_header_ics, "DATATYPE", "FLOAT"); - cb_data[i].guppiraw_header_ics.metadata.datashape.block_size = ctx.h_pwrbuf_size[i]/ctx.Nant; cb_data[i].guppiraw_header_ics.metadata.datashape.n_ant = 1; - cb_data[i].guppiraw_header.metadata.datashape.n_bit = (sizeof(float)*8)/2; // to account for assumption of data being complex + cb_data[i].guppiraw_header_ics.metadata.datashape.n_bit = sizeof(float)*8; cb_data[i].guppiraw_header_ics.metadata.datashape.n_time = ctx.Nds[i]; - cb_data[i].guppiraw_header_ics.metadata.datashape.n_pol = ctx.Npolout[i]; cb_data[i].guppiraw_header_ics.metadata.datashape.n_aspectchan = ctx.Nts[i] * ctx.Nc/ctx.Nant; guppiraw_header_put_metadata(&cb_data[i].guppiraw_header); } @@ -827,9 +845,8 @@ int main(int argc, char *argv[]) if(cb_data[i].per_ant_out) { cb_data[i].guppiraw_header.metadata.datashape.n_ant = 1; } - cb_data[i].guppiraw_header.metadata.datashape.n_bit = (sizeof(float)*8)/2; // to account for assumption of data being complex + cb_data[i].guppiraw_header.metadata.datashape.n_bit = sizeof(float)*8; cb_data[i].guppiraw_header.metadata.datashape.n_time = ctx.Nds[i]; - cb_data[i].guppiraw_header.metadata.datashape.n_pol = ctx.Npolout[i]; cb_data[i].guppiraw_header.metadata.datashape.n_aspectchan = ctx.Nts[i] * ctx.Nc/ctx.Nant; guppiraw_header_put_string(&cb_data[i].guppiraw_header, "DATATYPE", "FLOAT"); guppiraw_header_put_metadata(&cb_data[i].guppiraw_header); diff --git a/rawspec.h b/rawspec.h index e09c45f..ed73293 100644 --- a/rawspec.h +++ b/rawspec.h @@ -114,6 +114,9 @@ struct rawspec_context_s { // the next call to rawspec_initialize(). int input_conjugated; + // Flag indicating that the output is complex, a direct output of the FFT. + int complex_output; + // Flag indicating the concurrent output of the output data's incoherent-sum. int incoherently_sum; int Naws; diff --git a/rawspec_file.c b/rawspec_file.c index f9805fa..b142782 100644 --- a/rawspec_file.c +++ b/rawspec_file.c @@ -127,39 +127,44 @@ void * dump_file_thread_func(void *arg) { callback_data_t * cb_data = (callback_data_t *)arg; int retcode; - const size_t spectra_stride = cb_data->h_pwrbuf_size / (cb_data->Nds * sizeof(float)); - const size_t pol_stride = spectra_stride / cb_data->fb_hdr.nifs; - const size_t ant_stride = pol_stride / cb_data->Nant; - const size_t chan_stride = ant_stride / cb_data->Nf; // Single antenna case if(cb_data->fd && cb_data->h_pwrbuf) { if(cb_data->flag_file_output == FILE_FORMAT_GUPPIRAW) { + const char sample_bytesize = 2 * sizeof(float); + const size_t spectra_stride = cb_data->h_pwrbuf_size / (cb_data->Nds * sample_bytesize); + const size_t pol_stride = spectra_stride / cb_data->guppiraw_header.metadata.datashape.n_pol; + const size_t ant_stride = pol_stride / cb_data->Nant; + const size_t chan_stride = pol_stride / cb_data->Nf; // Nf is OBSNCHAN + int64_t bytes_written; + // printf("guppiraw_write antenna: %ld bytes\n\t", cb_data->h_pwrbuf_size / (cb_data->per_ant_out ? cb_data->Nant : 1 )); + fprintf(stdout, "guppiraw_write antenna start: %ld bytes\n", cb_data->h_pwrbuf_size); for(size_t i = 0; i < (cb_data->per_ant_out ? cb_data->Nant : 1); i++){ // printf("%d ", cb_data->fd[i]); if(cb_data->fd[i] == -1){ // Assume that the following file-descriptors aren't valid break; } - size_t bytes_written = guppiraw_write_block_arbitrary( + bytes_written = guppiraw_write_block_arbitrary( cb_data->fd[i], &cb_data->guppiraw_header, - cb_data->h_pwrbuf + ant_stride*i, - ant_stride*sizeof(float), - chan_stride*sizeof(float), - spectra_stride*sizeof(float), - pol_stride*sizeof(float) + cb_data->h_pwrbuf + ant_stride*2*i, + ant_stride*sample_bytesize, + chan_stride*sample_bytesize, + spectra_stride*sample_bytesize, + pol_stride*sample_bytesize ); if(bytes_written - < cb_data->h_pwrbuf_size / (cb_data->per_ant_out ? cb_data->Nant : 1) + < (int64_t) cb_data->h_pwrbuf_size / (cb_data->per_ant_out ? cb_data->Nant : 1) ) { cb_data->exit_soon = 1; cb_data->output_thread_valid = 0; fprintf(stderr, "GUPPIRAW-WRITE-ERROR: wrote %ld (%ld vs %ld)\n", bytes_written, cb_data->h_pwrbuf_size, cb_data->guppiraw_header.metadata.datashape.block_size); + break; } } - // printf("\n"); + fprintf(stdout, "guppiraw_write antenna end\n"); } else if(cb_data->Nant == 1) { if(cb_data->debug_callback) @@ -188,6 +193,9 @@ void * dump_file_thread_func(void *arg) break; } } else if(cb_data->per_ant_out && cb_data->Nant > 1) { + const size_t spectra_stride = cb_data->h_pwrbuf_size / (cb_data->Nds * sizeof(float)); + const size_t pol_stride = spectra_stride / cb_data->fb_hdr.nifs; + const size_t ant_stride = pol_stride / cb_data->Nant; // Multiple antennas, split output for(size_t k = 0; k < cb_data->Nds; k++){// Spectra out for(size_t j = 0; j < cb_data->fb_hdr.nifs; j++){// Npolout @@ -232,8 +240,7 @@ void * dump_file_thread_func(void *arg) if(cb_data->debug_callback) printf("dump_file_thread_func: write for ICS\n"); - switch(cb_data->flag_file_output) { - case FILE_FORMAT_FBH5: + if (cb_data->flag_file_output == FILE_FORMAT_FBH5) { retcode = fbh5_write(&(cb_data->fbh5_ctx_ics), &(cb_data->fb_hdr), cb_data->h_icsbuf, @@ -243,8 +250,7 @@ void * dump_file_thread_func(void *arg) cb_data->exit_soon = 1; cb_data->output_thread_valid = 0; } - break; - case FILE_FORMAT_FBSIGPROC: + } else if (cb_data->flag_file_output == FILE_FORMAT_FBSIGPROC) { if(write(cb_data->fd_ics, cb_data->h_icsbuf, cb_data->h_pwrbuf_size/cb_data->Nant) < 0) { @@ -252,24 +258,27 @@ void * dump_file_thread_func(void *arg) cb_data->output_thread_valid = 0; fprintf(stderr, "SIGPROC-WRITE-ERROR\n"); } - break; - case FILE_FORMAT_GUPPIRAW: + } else if (cb_data->flag_file_output == FILE_FORMAT_GUPPIRAW) { + const char sample_bytesize = 2 * sizeof(float); + const size_t spectra_stride = cb_data->h_pwrbuf_size / (cb_data->Nds * sample_bytesize); + const size_t pol_stride = spectra_stride / cb_data->guppiraw_header.metadata.datashape.n_pol; + const size_t ant_stride = pol_stride / cb_data->Nant; + const size_t chan_stride = pol_stride / cb_data->Nf; // Nf is OBSNCHAN if( guppiraw_write_block_arbitrary( cb_data->fd_ics, &cb_data->guppiraw_header_ics, cb_data->h_pwrbuf, - ant_stride*sizeof(float), - chan_stride*sizeof(float), - spectra_stride*sizeof(float), - pol_stride*sizeof(float) - ) < cb_data->h_pwrbuf_size / (cb_data->per_ant_out ? 1 : cb_data->Nant) + ant_stride*sample_bytesize, + chan_stride*sample_bytesize, + spectra_stride*sample_bytesize, + pol_stride*sample_bytesize + ) < (int64_t) cb_data->h_pwrbuf_size / (cb_data->per_ant_out ? 1 : cb_data->Nant) ) { cb_data->exit_soon = 1; cb_data->output_thread_valid = 0; fprintf(stderr, "GUPPIRAW-WRITE-ERROR\n"); } - break; } } diff --git a/rawspec_gpu.cu b/rawspec_gpu.cu index a5c2acc..1eb0581 100644 --- a/rawspec_gpu.cu +++ b/rawspec_gpu.cu @@ -160,6 +160,32 @@ __device__ void store_callback(void *p_v_out, ((float *)p_v_user)[offset] += pwr; } +// For complex-only mode (Npolout == 2), the store_callback_complex just needs to +// accumulate the elements into the power buffer. Pol0 is stored across the first half. +__device__ void store_callback_pol0_complex(void *p_v_out, + size_t offset, + cufftComplex element, + void *p_v_user, + void *p_v_shared) +{ + store_cb_data_t * d_scb_data = (store_cb_data_t *)p_v_user; + d_scb_data->pwr_buf_p00_i[2*offset] += element.x; + d_scb_data->pwr_buf_p00_i[2*offset + 1] += element.y; +} + +// For complex-only mode (Npolout == 2), the store_callback_complex just needs to +// accumulate the elements into the power buffer. Pol1 is stored across the second half. +__device__ void store_callback_pol1_complex(void *p_v_out, + size_t offset, + cufftComplex element, + void *p_v_user, + void *p_v_shared) +{ + store_cb_data_t * d_scb_data = (store_cb_data_t *)p_v_user; + d_scb_data->pwr_buf_p01_re_u[2*offset] += element.x; + d_scb_data->pwr_buf_p01_re_u[2*offset + 1] += element.y; +} + // For full-Stokes mode, the store_callback_pol0_iquv function stores the // voltage data into the first half of the 2x-sized FFT output buffer and // accummulates (i.e. adds) the pol0 power into the first two quarters (I and @@ -268,6 +294,8 @@ __device__ void store_callback_pol1_conj(void *p_v_out, __device__ cufftCallbackLoadC d_cufft_load_callback = load_callback; __device__ cufftCallbackStoreC d_cufft_store_callback = store_callback; +__device__ cufftCallbackStoreC d_cufft_store_callback_complex_pol0 = store_callback_pol0_complex; +__device__ cufftCallbackStoreC d_cufft_store_callback_complex_pol1 = store_callback_pol1_complex; __device__ cufftCallbackStoreC d_cufft_store_callback_pol0 = store_callback_pol0; __device__ cufftCallbackStoreC d_cufft_store_callback_pol1 = store_callback_pol1; __device__ cufftCallbackStoreC d_cufft_store_callback_pol1_conj = store_callback_pol1_conj; @@ -284,12 +312,12 @@ __device__ cufftCallbackStoreC d_cufft_store_callback_pol1_iquv_conj = store_cal // Expectation of blockDim, with 1 thread each: // grid.x = width; // grid.y = height; -// grid.z = 1; +// grid.z = complexity_factor; // All pitches are in units of float __global__ void mem2dAccumMoveFloat(float* dst, size_t dpitch, float* src, size_t spitch) { - float* src_offst = src + blockDim.x*blockIdx.y*spitch + blockIdx.x; - dst[blockDim.x*blockIdx.y*dpitch + blockIdx.x] += *src_offst; + float* src_offst = src + blockIdx.z + blockDim.z*(blockDim.x*blockIdx.y*spitch + blockIdx.x); + dst[blockIdx.z + blockDim.z*(blockDim.x*blockIdx.y*dpitch + blockIdx.x)] += *src_offst; *src_offst = 0.0; } @@ -465,6 +493,7 @@ int rawspec_initialize(rawspec_context * ctx) // Host copies of cufft callback pointers cufftCallbackLoadC h_cufft_load_callback; cufftCallbackStoreC h_cufft_store_callback; + cufftCallbackStoreC h_cufft_store_callback_complex_pols[2]; cufftCallbackStoreC h_cufft_store_callback_pols[2]; cufftCallbackStoreC h_cufft_store_callback_iquv[2]; @@ -486,8 +515,22 @@ int rawspec_initialize(rawspec_context * ctx) // Validate/set Npolout values for(i=0; iNo; i++) { - if(abs(ctx->Npolout[i]) != 4 || ctx->Np != 2) { - ctx->Npolout[i] = 1; + switch(ctx->Npolout[i]) { + case -4: + case 4: + if(ctx->Np != 2) { + ctx->Npolout[i] = 1; + } + break; + + case -2: + case -1: + ctx->complex_output = 1; + ctx->Npolout[i] = abs(ctx->Npolout[i]); + break; + + default: + ctx->Npolout[i] = 1; } } @@ -747,17 +790,20 @@ int rawspec_initialize(rawspec_context * ctx) } // Calculate grid dimensions - gpu_ctx->grid[i].x = (ctx->Nts[i] + MAX_THREADS - 1) / MAX_THREADS; + gpu_ctx->nthreads[i] = ctx->Nts[i]; + gpu_ctx->nthreads[i] *= (ctx->complex_output + 1)%3; + gpu_ctx->grid[i].x = (gpu_ctx->nthreads[i] + MAX_THREADS - 1) / MAX_THREADS; gpu_ctx->grid[i].y = ctx->Nds[i]; gpu_ctx->grid[i].z = ctx->Nbc; // Calculate number of threads per block - gpu_ctx->nthreads[i] = MIN(ctx->Nts[i], MAX_THREADS); + gpu_ctx->nthreads[i] = MIN(gpu_ctx->nthreads[i], MAX_THREADS); // Host buffer needs to accommodate the number of integrations that will be // dumped at one time (Nd). ctx->h_pwrbuf_size[i] = abs(ctx->Npolout[i]) * ctx->Nds[i]*ctx->Nts[i]*ctx->Nc*sizeof(float); + ctx->h_pwrbuf_size[i] *= (1 + ctx->complex_output)%3; // complex-output conditional double #ifdef VERBOSE_ALLOC printf("FFT Host dump buffer[%d] size == %lu\n", i, ctx->h_pwrbuf_size[i]); #endif @@ -949,13 +995,15 @@ int rawspec_initialize(rawspec_context * ctx) // For each output product for(i=0; i < ctx->No; i++) { // Power output buffer + buf_size = abs(ctx->Npolout[i]) * ctx->Nb*ctx->Ntpb*ctx->Nbc*sizeof(float); + buf_size *= (1 + ctx->complex_output)%3; // complex-output conditional double #ifdef VERBOSE_ALLOC - printf("Power output buffer size == %u * %lu == %lu\n", - abs(ctx->Npolout[i]), ctx->Nb*ctx->Ntpb*ctx->Nbc*sizeof(float), - abs(ctx->Npolout[i]) * ctx->Nb*ctx->Ntpb*ctx->Nbc*sizeof(float)); + printf("Power output buffer size == abs(%u) * %lu == %lu\n", + ctx->Npolout[i], ctx->Nb*ctx->Ntpb*ctx->Nbc*sizeof(float), + buf_size); #endif cuda_rc = cudaMalloc(&gpu_ctx->d_pwr_out[i], - abs(ctx->Npolout[i]) * ctx->Nb*ctx->Ntpb*ctx->Nbc*sizeof(float)); + buf_size); if(cuda_rc != cudaSuccess) { PRINT_CUDA_ERRMSG(cuda_rc); rawspec_cleanup(ctx); @@ -963,7 +1011,7 @@ int rawspec_initialize(rawspec_context * ctx) } // Clear power output buffer cuda_rc = cudaMemset(gpu_ctx->d_pwr_out[i], 0, - abs(ctx->Npolout[i]) * ctx->Nb*ctx->Ntpb*ctx->Nbc*sizeof(float)); + buf_size); if(cuda_rc != cudaSuccess) { PRINT_CUDA_ERRMSG(cuda_rc); rawspec_cleanup(ctx); @@ -971,13 +1019,15 @@ int rawspec_initialize(rawspec_context * ctx) } if(gpu_ctx->Nis[i] > 1 && ctx->Nbc < ctx->Nc){ // Full power output buffer + buf_size = abs(ctx->Npolout[i]) * ctx->Nb*ctx->Ntpb*ctx->Nc*sizeof(float); + buf_size *= (1 + ctx->complex_output)%3; // complex-output conditional double #ifdef VERBOSE_ALLOC printf("Full power output buffer cache size == %u * %lu == %lu\n", - abs(ctx->Npolout[i]), ctx->Nb*ctx->Ntpb*ctx->Nc*sizeof(float), - abs(ctx->Npolout[i]) * ctx->Nb*ctx->Ntpb*ctx->Nc*sizeof(float)); + ctx->Npolout[i], ctx->Nb*ctx->Ntpb*ctx->Nc*sizeof(float), + buf_size); #endif cuda_rc = cudaMalloc(&gpu_ctx->d_prev_pwr_out_cache[i], - abs(ctx->Npolout[i]) * ctx->Nb*ctx->Ntpb*ctx->Nc*sizeof(float)); + buf_size); if(cuda_rc != cudaSuccess) { PRINT_CUDA_ERRMSG(cuda_rc); rawspec_cleanup(ctx); @@ -985,7 +1035,7 @@ int rawspec_initialize(rawspec_context * ctx) } // Clear power output buffer cuda_rc = cudaMemset(gpu_ctx->d_prev_pwr_out_cache[i], 0, - abs(ctx->Npolout[i]) * ctx->Nb*ctx->Ntpb*ctx->Nc*sizeof(float)); + buf_size); if(cuda_rc != cudaSuccess) { PRINT_CUDA_ERRMSG(cuda_rc); rawspec_cleanup(ctx); @@ -997,13 +1047,15 @@ int rawspec_initialize(rawspec_context * ctx) } if(ctx->incoherently_sum){ + buf_size = abs(ctx->Npolout[i]) * ctx->Nb*ctx->Ntpb*ctx->Nc*sizeof(float)/ctx->Nant; + buf_size *= (1 + ctx->complex_output)%3; // complex-output conditional double #ifdef VERBOSE_ALLOC printf("ICS output buffer size == %u * %lu / %u == %lu\n", - abs(ctx->Npolout[i]), ctx->Nb*ctx->Ntpb*ctx->Nc*sizeof(float), ctx->Nant, - abs(ctx->Npolout[i]) * ctx->Nb*ctx->Ntpb*ctx->Nc*sizeof(float)/ctx->Nant); + ctx->Npolout[i], ctx->Nb*ctx->Ntpb*ctx->Nc*sizeof(float), ctx->Nant, + buf_size); #endif cuda_rc = cudaMalloc(&gpu_ctx->d_ics_out[i], - abs(ctx->Npolout[i]) * ctx->Nb*ctx->Ntpb*ctx->Nc*sizeof(float)/ctx->Nant); + buf_size); if(cuda_rc != cudaSuccess) { PRINT_CUDA_ERRMSG(cuda_rc); rawspec_cleanup(ctx); @@ -1011,7 +1063,7 @@ int rawspec_initialize(rawspec_context * ctx) } // Clear incoherent-sum output buffer cuda_rc = cudaMemset(gpu_ctx->d_ics_out[i], 0, - abs(ctx->Npolout[i]) * ctx->Nb*ctx->Ntpb*ctx->Nc*sizeof(float)/ctx->Nant); + buf_size); if(cuda_rc != cudaSuccess) { PRINT_CUDA_ERRMSG(cuda_rc); rawspec_cleanup(ctx); @@ -1100,6 +1152,24 @@ int rawspec_initialize(rawspec_context * ctx) return 1; } + cuda_rc = cudaMemcpyFromSymbol(&h_cufft_store_callback_complex_pols[0], + d_cufft_store_callback_complex_pol0, + sizeof(h_cufft_store_callback_complex_pols[0])); + if(cuda_rc != cudaSuccess) { + PRINT_CUDA_ERRMSG(cuda_rc); + rawspec_cleanup(ctx); + return 1; + } + + cuda_rc = cudaMemcpyFromSymbol(&h_cufft_store_callback_complex_pols[1], + d_cufft_store_callback_complex_pol1, + sizeof(h_cufft_store_callback_complex_pols[1])); + if(cuda_rc != cudaSuccess) { + PRINT_CUDA_ERRMSG(cuda_rc); + rawspec_cleanup(ctx); + return 1; + } + cuda_rc = cudaMemcpyFromSymbol(&h_cufft_store_callback_pols[0], d_cufft_store_callback_pol0, sizeof(h_cufft_store_callback_pols[0])); @@ -1207,7 +1277,12 @@ int rawspec_initialize(rawspec_context * ctx) return 1; } // Store callback(s) - if(ctx->Npolout[i] == 1) { + if(ctx->complex_output == 1) { + cufft_rc = cufftXtSetCallback(gpu_ctx->plan[i][p], + (void **)&h_cufft_store_callback_complex_pols[p], + CUFFT_CB_ST_COMPLEX, + (void **)&gpu_ctx->d_scb_data[i]); + } else if(ctx->Npolout[i] == 1) { cufft_rc = cufftXtSetCallback(gpu_ctx->plan[i][p], (void **)&h_cufft_store_callback, CUFFT_CB_ST_COMPLEX, @@ -1525,12 +1600,14 @@ int rawspec_start_processing(rawspec_context * ctx, int fft_dir) size_t fft_outbuf_length; dim3 grid_ics; dim3 grid_full_pwr; - grid_full_pwr.x = ctx->Nb*ctx->Ntpb*ctx->Nbc; - grid_full_pwr.z = 1; const size_t Nchan_per_antenna = ctx->Nc/ctx->Nant; const size_t Nantenna_per_batch = ctx->Nbc/Nchan_per_antenna; + const char complexity_factor = ctx->complex_output ? 2 : 1; bool is_last_channel_batch; + grid_full_pwr.x = ctx->Nb*ctx->Ntpb*ctx->Nbc; + grid_full_pwr.z = complexity_factor; + // Increment inbuf_count gpu_ctx->inbuf_count++; @@ -1568,9 +1645,9 @@ int rawspec_start_processing(rawspec_context * ctx, int fft_dir) // mem2dAccumMoveFloat memset(0)s the source (d_prev_pwr_out_cache) mem2dAccumMoveFloat<<compute_stream >>>(gpu_ctx->d_pwr_out[i], - ctx->Nb*ctx->Ntpb*ctx->Nbc, - gpu_ctx->d_prev_pwr_out_cache[i] + (c * abs(ctx->Npolout[i]) * ctx->Nb*ctx->Ntpb), - ctx->Nb*ctx->Ntpb*ctx->Nc + ctx->Nb*ctx->Ntpb*ctx->Nbc*complexity_factor, + gpu_ctx->d_prev_pwr_out_cache[i] + (c * abs(ctx->Npolout[i]) * ctx->Nb*ctx->Ntpb * complexity_factor), + ctx->Nb*ctx->Ntpb*ctx->Nc*complexity_factor ); } @@ -1583,11 +1660,11 @@ int rawspec_start_processing(rawspec_context * ctx, int fft_dir) gpu_ctx->nthreads[i], 0, gpu_ctx->compute_stream>>> ( - gpu_ctx->d_pwr_out[i] + p*ctx->Nb*ctx->Ntpb*ctx->Nbc, + gpu_ctx->d_pwr_out[i] + p*ctx->Nb*ctx->Ntpb*ctx->Nbc*complexity_factor, MIN(ctx->Nas[i], gpu_ctx->Nss[i]), // Na - ctx->Nts[i], // xpitch - ctx->Nas[i]*ctx->Nts[i], // ypitch - ctx->Nb*ctx->Ntpb // zpitch + complexity_factor*ctx->Nts[i], // xpitch + complexity_factor*ctx->Nas[i]*ctx->Nts[i], // ypitch + complexity_factor*ctx->Nb*ctx->Ntpb // zpitch ); } } @@ -1640,16 +1717,16 @@ int rawspec_start_processing(rawspec_context * ctx, int fft_dir) // Copy integrated power spectra (or spectrum) to host. This is done as // two 2D copies to get channel 0 in the center of the spectrum. Special // care is taken in the unlikely event that Nt is odd. - src = gpu_ctx->d_pwr_out[i] + p*ctx->Nb*ctx->Ntpb*ctx->Nbc; - dst = ctx->h_pwrbuf[i] + (p*ctx->Nts[i]*ctx->Nc) + (c*ctx->Nts[i]); - spitch = gpu_ctx->Nss[i] * ctx->Nts[i] * sizeof(float); - dpitch = ctx->Nts[i] * sizeof(float); + src = gpu_ctx->d_pwr_out[i] + (p*ctx->Nb*ctx->Ntpb*ctx->Nbc)*complexity_factor; + dst = ctx->h_pwrbuf[i] + ((p*ctx->Nts[i]*ctx->Nc) + (c*ctx->Nts[i]))*complexity_factor; + dpitch = ctx->Nts[i] * sizeof(float) * complexity_factor; + spitch = gpu_ctx->Nss[i] * dpitch; height = ctx->Nbc; for(d=0; d < ctx->Nds[i]; d++) { // Lo to hi - width = ((ctx->Nts[i]+1) / 2) * sizeof(float); + width = ((ctx->Nts[i]+1) / 2) * sizeof(float) * complexity_factor; cuda_rc = cudaMemcpy2DAsync(dst + ctx->Nts[i]/2, dpitch, src, @@ -1666,7 +1743,7 @@ int rawspec_start_processing(rawspec_context * ctx, int fft_dir) } // Hi to lo - width = (ctx->Nts[i] / 2) * sizeof(float); + width = (ctx->Nts[i] / 2) * sizeof(float) * complexity_factor; cuda_rc = cudaMemcpy2DAsync(dst, dpitch, src + (ctx->Nts[i]+1) / 2, From eccdbb415056896ef0049124a48bd62372ff6833 Mon Sep 17 00:00:00 2001 From: radonnachie Date: Wed, 3 Aug 2022 12:24:06 +0000 Subject: [PATCH 6/9] * complex-factor conditional doubling --- rawspec_gpu.cu | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/rawspec_gpu.cu b/rawspec_gpu.cu index 1eb0581..772496e 100644 --- a/rawspec_gpu.cu +++ b/rawspec_gpu.cu @@ -533,6 +533,7 @@ int rawspec_initialize(rawspec_context * ctx) ctx->Npolout[i] = 1; } } + ctx->complex_output = min(ctx->complex_output, 2); // Validate Ntpb if(ctx->Ntpb == 0) { @@ -791,7 +792,7 @@ int rawspec_initialize(rawspec_context * ctx) // Calculate grid dimensions gpu_ctx->nthreads[i] = ctx->Nts[i]; - gpu_ctx->nthreads[i] *= (ctx->complex_output + 1)%3; + gpu_ctx->nthreads[i] *= 1 + ctx->complex_output; gpu_ctx->grid[i].x = (gpu_ctx->nthreads[i] + MAX_THREADS - 1) / MAX_THREADS; gpu_ctx->grid[i].y = ctx->Nds[i]; gpu_ctx->grid[i].z = ctx->Nbc; @@ -803,7 +804,7 @@ int rawspec_initialize(rawspec_context * ctx) // dumped at one time (Nd). ctx->h_pwrbuf_size[i] = abs(ctx->Npolout[i]) * ctx->Nds[i]*ctx->Nts[i]*ctx->Nc*sizeof(float); - ctx->h_pwrbuf_size[i] *= (1 + ctx->complex_output)%3; // complex-output conditional double + ctx->h_pwrbuf_size[i] *= 1 + ctx->complex_output; // complex-output conditional double #ifdef VERBOSE_ALLOC printf("FFT Host dump buffer[%d] size == %lu\n", i, ctx->h_pwrbuf_size[i]); #endif @@ -996,7 +997,7 @@ int rawspec_initialize(rawspec_context * ctx) for(i=0; i < ctx->No; i++) { // Power output buffer buf_size = abs(ctx->Npolout[i]) * ctx->Nb*ctx->Ntpb*ctx->Nbc*sizeof(float); - buf_size *= (1 + ctx->complex_output)%3; // complex-output conditional double + buf_size *= 1 + ctx->complex_output; // complex-output conditional double #ifdef VERBOSE_ALLOC printf("Power output buffer size == abs(%u) * %lu == %lu\n", ctx->Npolout[i], ctx->Nb*ctx->Ntpb*ctx->Nbc*sizeof(float), @@ -1020,7 +1021,7 @@ int rawspec_initialize(rawspec_context * ctx) if(gpu_ctx->Nis[i] > 1 && ctx->Nbc < ctx->Nc){ // Full power output buffer buf_size = abs(ctx->Npolout[i]) * ctx->Nb*ctx->Ntpb*ctx->Nc*sizeof(float); - buf_size *= (1 + ctx->complex_output)%3; // complex-output conditional double + buf_size *= 1 + ctx->complex_output; // complex-output conditional double #ifdef VERBOSE_ALLOC printf("Full power output buffer cache size == %u * %lu == %lu\n", ctx->Npolout[i], ctx->Nb*ctx->Ntpb*ctx->Nc*sizeof(float), @@ -1048,7 +1049,7 @@ int rawspec_initialize(rawspec_context * ctx) if(ctx->incoherently_sum){ buf_size = abs(ctx->Npolout[i]) * ctx->Nb*ctx->Ntpb*ctx->Nc*sizeof(float)/ctx->Nant; - buf_size *= (1 + ctx->complex_output)%3; // complex-output conditional double + buf_size *= 1 + ctx->complex_output; // complex-output conditional double #ifdef VERBOSE_ALLOC printf("ICS output buffer size == %u * %lu / %u == %lu\n", ctx->Npolout[i], ctx->Nb*ctx->Ntpb*ctx->Nc*sizeof(float), ctx->Nant, @@ -1602,7 +1603,7 @@ int rawspec_start_processing(rawspec_context * ctx, int fft_dir) dim3 grid_full_pwr; const size_t Nchan_per_antenna = ctx->Nc/ctx->Nant; const size_t Nantenna_per_batch = ctx->Nbc/Nchan_per_antenna; - const char complexity_factor = ctx->complex_output ? 2 : 1; + const char complexity_factor = 1 + ctx->complex_output; bool is_last_channel_batch; grid_full_pwr.x = ctx->Nb*ctx->Ntpb*ctx->Nbc; From 293be667ddd077a785726085adbf07ea337409e7 Mon Sep 17 00:00:00 2001 From: radonnachie Date: Wed, 3 Aug 2022 12:25:03 +0000 Subject: [PATCH 7/9] * h_scb_data pointers use complex-factor --- rawspec_gpu.cu | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/rawspec_gpu.cu b/rawspec_gpu.cu index 772496e..25350e2 100644 --- a/rawspec_gpu.cu +++ b/rawspec_gpu.cu @@ -182,8 +182,8 @@ __device__ void store_callback_pol1_complex(void *p_v_out, void *p_v_shared) { store_cb_data_t * d_scb_data = (store_cb_data_t *)p_v_user; - d_scb_data->pwr_buf_p01_re_u[2*offset] += element.x; - d_scb_data->pwr_buf_p01_re_u[2*offset + 1] += element.y; + d_scb_data->pwr_buf_p11_q[2*offset] += element.x; + d_scb_data->pwr_buf_p11_q[2*offset + 1] += element.y; } // For full-Stokes mode, the store_callback_pol0_iquv function stores the @@ -1107,12 +1107,14 @@ int rawspec_initialize(rawspec_context * ctx) // so we can initialize them that way even if Npolout == 1 // (because they will never be used). It might be slightly // safer to init them to the same as pwr_buf_p00_i if Npolout == 1. + buf_size = ctx->Nb*ctx->Ntpb*ctx->Nbc; + buf_size *= 1 + ctx->complex_output; // complex-output conditional double h_scb_data.pwr_buf_p11_q = - gpu_ctx->d_pwr_out[i] + 1*ctx->Nb*ctx->Ntpb*ctx->Nbc; + gpu_ctx->d_pwr_out[i] + 1*buf_size; h_scb_data.pwr_buf_p01_re_u = - gpu_ctx->d_pwr_out[i] + 2*ctx->Nb*ctx->Ntpb*ctx->Nbc; + gpu_ctx->d_pwr_out[i] + 2*buf_size; h_scb_data.pwr_buf_p01_im_v = - gpu_ctx->d_pwr_out[i] + 3*ctx->Nb*ctx->Ntpb*ctx->Nbc; + gpu_ctx->d_pwr_out[i] + 3*buf_size; // Allocate device memory for store_cb_data_t array cuda_rc = cudaMalloc(&gpu_ctx->d_scb_data[i], sizeof(store_cb_data_t)); From 939d54c09cd7166c20ed046b1fde679b18cfd315 Mon Sep 17 00:00:00 2001 From: radonnachie Date: Wed, 3 Aug 2022 12:25:47 +0000 Subject: [PATCH 8/9] ^ GUPPI RAW output (at least for -t1 -f1) --- rawspec_gpu.cu | 88 +++++++++++++++++++++++++++++++------------------- 1 file changed, 54 insertions(+), 34 deletions(-) diff --git a/rawspec_gpu.cu b/rawspec_gpu.cu index 25350e2..182a28d 100644 --- a/rawspec_gpu.cu +++ b/rawspec_gpu.cu @@ -1727,44 +1727,64 @@ int rawspec_start_processing(rawspec_context * ctx, int fft_dir) height = ctx->Nbc; for(d=0; d < ctx->Nds[i]; d++) { - - // Lo to hi - width = ((ctx->Nts[i]+1) / 2) * sizeof(float) * complexity_factor; - cuda_rc = cudaMemcpy2DAsync(dst + ctx->Nts[i]/2, - dpitch, - src, - spitch, - width, - height, - cudaMemcpyDeviceToHost, - gpu_ctx->compute_stream); - - if(cuda_rc != cudaSuccess) { - PRINT_CUDA_ERRMSG(cuda_rc); - rawspec_cleanup(ctx); - return 1; + + if(ctx->complex_output) { + // GUPPI RAW output, don't translate channels + width = dpitch; + cuda_rc = cudaMemcpy2DAsync(dst, + dpitch, + src, + spitch, + width, + height, + cudaMemcpyDeviceToHost, + gpu_ctx->compute_stream); + + if(cuda_rc != cudaSuccess) { + PRINT_CUDA_ERRMSG(cuda_rc); + rawspec_cleanup(ctx); + return 1; + } } - - // Hi to lo - width = (ctx->Nts[i] / 2) * sizeof(float) * complexity_factor; - cuda_rc = cudaMemcpy2DAsync(dst, - dpitch, - src + (ctx->Nts[i]+1) / 2, - spitch, - width, - height, - cudaMemcpyDeviceToHost, - gpu_ctx->compute_stream); - - if(cuda_rc != cudaSuccess) { - PRINT_CUDA_ERRMSG(cuda_rc); - rawspec_cleanup(ctx); - return 1; + else { + // Lo to hi + width = ((ctx->Nts[i]+1) / 2) * sizeof(float); + cuda_rc = cudaMemcpy2DAsync(dst + ctx->Nts[i]/2, + dpitch, + src, + spitch, + width, + height, + cudaMemcpyDeviceToHost, + gpu_ctx->compute_stream); + + if(cuda_rc != cudaSuccess) { + PRINT_CUDA_ERRMSG(cuda_rc); + rawspec_cleanup(ctx); + return 1; + } + + // Hi to lo + width = (ctx->Nts[i] / 2) * sizeof(float); + cuda_rc = cudaMemcpy2DAsync(dst, + dpitch, + src + (ctx->Nts[i]+1) / 2, + spitch, + width, + height, + cudaMemcpyDeviceToHost, + gpu_ctx->compute_stream); + + if(cuda_rc != cudaSuccess) { + PRINT_CUDA_ERRMSG(cuda_rc); + rawspec_cleanup(ctx); + return 1; + } } // Increment src and dst pointers - src += ctx->Nts[i] * ctx->Nas[i]; - dst += abs(ctx->Npolout[i]) * ctx->Nts[i] * ctx->Nc; + src += ctx->Nts[i] * ctx->Nas[i] * complexity_factor; + dst += abs(ctx->Npolout[i]) * ctx->Nts[i] * ctx->Nc * complexity_factor; } } From 6aadb76ff6df8c1564d52d862273142674e6c360 Mon Sep 17 00:00:00 2001 From: radonnachie Date: Wed, 24 Aug 2022 12:35:02 +0000 Subject: [PATCH 9/9] !! WIP --- compare.jl | 77 +++++++++++++++++++++++++++++++ plot.jl | 36 +++++++++++++++ rawspec.c | 5 +- rawspec_file.c | 44 ++++-------------- rawspec_gpu.cu | 112 ++++++++++++++++++++++++++++++++++----------- rawspec_rawutils.c | 2 +- 6 files changed, 213 insertions(+), 63 deletions(-) create mode 100644 compare.jl create mode 100644 plot.jl diff --git a/compare.jl b/compare.jl new file mode 100644 index 0000000..76bc53d --- /dev/null +++ b/compare.jl @@ -0,0 +1,77 @@ +using Printf +using FFTW +using Blio: GuppiRaw + +struct Result + message::Union{String, Nothing} + value::Bool +end + +function _fft(guppidata, points) + data = reshape(guppidata, (size(guppidata, 1), points, :, size(guppidata)[3:end]...)) + # [pol, fine_spec, time, coarse_chan, antenna] + spectra = fft(data, 2) + # [pol, fine_chan, time, antenna] + spectra = cat((spectra[:, :, :, i, :] for i in 1:size(spectra,4))..., dims=2) + # [pol, time, fine_chan, antenna] + return permutedims(spectra, [1,3,2,4]) +end + +function accumulate(guppidata, samples) + data = reshape(guppidata, (size(guppidata, 1), samples, :, size(guppidata)[3:end]...)) + acc = sum(data, dims=2) + return reshape(acc, (size(guppidata, 1), :, size(guppidata)[3:end]...)) +end + +function mapToFloat(value::Integer, type::Type) + return value < 0 ? -1.0*(value / typemin(type)) : value / typemax(type) +end + +function mapToFloat(value::Complex{<:Integer}, type::Type) + return complex(mapToFloat(real(value), real(type)), mapToFloat(imag(value), real(type))) +end + +function compare(i_data, o_data, atol=0.01)::Result + if size(i_data) != size(o_data) + return Result(@sprintf("Shape mismatch: %s != %s", size(i_data), size(o_data)), false) + end + dims_correct = Array{Bool}(undef, size(i_data)[2:end]) + for i in CartesianIndices(dims_correct) + dims_correct[i] = all(isapprox.(real(i_data[:, i]), real(o_data[:, i]), atol=atol)) && all(isapprox.(imag(i_data[:, i]), imag(o_data[:, i]), atol=atol)) + if !dims_correct[i] + println(Result(@sprintf("Pol data mismatch @ %s: %s != %s\n\t(atol: %s)", i, i_data[:, i], o_data[:, i], i_data[:, i] - o_data[:, i]), false)) + end + end + + Result(nothing, all(dims_correct)) +end + +i_grheader = GuppiRaw.Header() +o_grheader = GuppiRaw.Header() + +i_fio = open("/mnt/buf0/test.0000.raw", "r") +o_fio = open("/mnt/buf0/test.rawspec.0000.raw", "r") + + read!(i_fio, i_grheader) + i_data = Array(i_grheader) + read!(i_fio, i_data) + i_data = mapToFloat.(i_data, eltype(i_data)) + + read!(o_fio, o_grheader) + o_data = Array(o_grheader) + read!(o_fio, o_data) + + fftpoints = div(size(o_data, 3), size(i_data, 3)) + accumulation = div(div(size(i_data, 2), size(o_data, 2)), fftpoints) + println("fftpoints: ", fftpoints) + println("accumulation: ", accumulation) + + upchannelized = fftpoints == 1 ? i_data : _fft(i_data, fftpoints) + accum_fine = accumulation == 1 ? upchannelized : accumulate(upchannelized, accumulation) + + atol = 0.015 * max(log2(fftpoints), fftpoints) * accumulation + + println(compare(accum_fine, o_data, atol)) + +close(i_fio) +close(o_fio) \ No newline at end of file diff --git a/plot.jl b/plot.jl new file mode 100644 index 0000000..47a3c4e --- /dev/null +++ b/plot.jl @@ -0,0 +1,36 @@ +using Printf +using Plots +using Blio: GuppiRaw + +ENV["GKSwstype"]="nul" # disable plots display + +dir = "/mnt/buf0/" + +for filestem in ["bladetest_signal_out"] + + grheader = GuppiRaw.Header() + fio = open(@sprintf("%s%s.0000.raw", dir, filestem), "r") + read!(fio, grheader) + data = Array(grheader) + @printf("%s, %s bytes\n", size(grheader), prod(size(grheader))*2) + @printf("%d entries, %s bytes\n", length(grheader), length(grheader)*80) + @printf("%s\n", grheader) + @printf("@%d, %d bytes in file (%d bytes remain)\n", position(fio), filesize(fio), filesize(fio) - position(fio)) + read!(fio, data) + # read!(fio, grheader) + # read!(fio, data) + data = abs.(data) + @printf("data [%s, %s]\n", min(data...), max(data...)) + if min(data...) != NaN + + figpol1 = heatmap(data[1, :, 1:8, 1]) + figpol2 = heatmap(data[2, :, 1:8, 1]) + # title!(fig, "New title") + # flipxaxis!(fig, "New xlabel") + # yaxis!(fig, "New ylabel") + + savefig(plot(figpol1, figpol2, layout=(2,1)), @sprintf("%s.png", filestem)) + end + + close(fio) +end \ No newline at end of file diff --git a/rawspec.c b/rawspec.c index 4b42daf..878c158 100644 --- a/rawspec.c +++ b/rawspec.c @@ -829,7 +829,7 @@ int main(int argc, char *argv[]) guppiraw_header_copy(&cb_data[i].guppiraw_header_ics, &guppiraw_header); cb_data[i].guppiraw_header_ics.metadata.user_data = malloc(sizeof(rawspec_raw_hdr_t)); memcpy(cb_data[i].guppiraw_header_ics.metadata.user_data, &raw_hdr, sizeof(rawspec_raw_hdr_t)); - guppiraw_header_put_string(&cb_data[i].guppiraw_header_ics, "DATATYPE", "FLOAT"); + guppiraw_header_put_string(&cb_data[i].guppiraw_header_ics, "DATATYPE", "INTEGER"); cb_data[i].guppiraw_header_ics.metadata.datashape.n_ant = 1; cb_data[i].guppiraw_header_ics.metadata.datashape.n_bit = sizeof(float)*8; @@ -848,7 +848,8 @@ int main(int argc, char *argv[]) cb_data[i].guppiraw_header.metadata.datashape.n_bit = sizeof(float)*8; cb_data[i].guppiraw_header.metadata.datashape.n_time = ctx.Nds[i]; cb_data[i].guppiraw_header.metadata.datashape.n_aspectchan = ctx.Nts[i] * ctx.Nc/ctx.Nant; - guppiraw_header_put_string(&cb_data[i].guppiraw_header, "DATATYPE", "FLOAT"); + guppiraw_header_put_string(&cb_data[i].guppiraw_header, "DATATYPE", "INTEGER"); + guppiraw_header_put_double(&cb_data[i].guppiraw_header, "CHAN_BW", raw_hdr.obsbw/cb_data[i].guppiraw_header.metadata.datashape.n_aspectchan); guppiraw_header_put_metadata(&cb_data[i].guppiraw_header); char* header_string = guppiraw_header_malloc_string(&guppiraw_header); printf("```%s```\n", header_string); diff --git a/rawspec_file.c b/rawspec_file.c index b142782..6eec236 100644 --- a/rawspec_file.c +++ b/rawspec_file.c @@ -131,40 +131,25 @@ void * dump_file_thread_func(void *arg) // Single antenna case if(cb_data->fd && cb_data->h_pwrbuf) { if(cb_data->flag_file_output == FILE_FORMAT_GUPPIRAW) { - const char sample_bytesize = 2 * sizeof(float); - const size_t spectra_stride = cb_data->h_pwrbuf_size / (cb_data->Nds * sample_bytesize); - const size_t pol_stride = spectra_stride / cb_data->guppiraw_header.metadata.datashape.n_pol; - const size_t ant_stride = pol_stride / cb_data->Nant; - const size_t chan_stride = pol_stride / cb_data->Nf; // Nf is OBSNCHAN - int64_t bytes_written; + const size_t ant_stride = cb_data->h_pwrbuf_size / cb_data->Nant; - // printf("guppiraw_write antenna: %ld bytes\n\t", cb_data->h_pwrbuf_size / (cb_data->per_ant_out ? cb_data->Nant : 1 )); - fprintf(stdout, "guppiraw_write antenna start: %ld bytes\n", cb_data->h_pwrbuf_size); for(size_t i = 0; i < (cb_data->per_ant_out ? cb_data->Nant : 1); i++){ - // printf("%d ", cb_data->fd[i]); if(cb_data->fd[i] == -1){ // Assume that the following file-descriptors aren't valid break; } - bytes_written = guppiraw_write_block_arbitrary( - cb_data->fd[i], - &cb_data->guppiraw_header, - cb_data->h_pwrbuf + ant_stride*2*i, - ant_stride*sample_bytesize, - chan_stride*sample_bytesize, - spectra_stride*sample_bytesize, - pol_stride*sample_bytesize - ); - if(bytes_written - < (int64_t) cb_data->h_pwrbuf_size / (cb_data->per_ant_out ? cb_data->Nant : 1) + if(guppiraw_write_block( + cb_data->fd[i], + &cb_data->guppiraw_header, + cb_data->h_pwrbuf + ant_stride*i + ) < (int64_t) cb_data->h_pwrbuf_size / (cb_data->per_ant_out ? cb_data->Nant : 1) ) { cb_data->exit_soon = 1; cb_data->output_thread_valid = 0; - fprintf(stderr, "GUPPIRAW-WRITE-ERROR: wrote %ld (%ld vs %ld)\n", bytes_written, cb_data->h_pwrbuf_size, cb_data->guppiraw_header.metadata.datashape.block_size); + fprintf(stderr, "GUPPIRAW-WRITE-ERROR\n"); break; } } - fprintf(stdout, "guppiraw_write antenna end\n"); } else if(cb_data->Nant == 1) { if(cb_data->debug_callback) @@ -259,21 +244,12 @@ void * dump_file_thread_func(void *arg) fprintf(stderr, "SIGPROC-WRITE-ERROR\n"); } } else if (cb_data->flag_file_output == FILE_FORMAT_GUPPIRAW) { - const char sample_bytesize = 2 * sizeof(float); - const size_t spectra_stride = cb_data->h_pwrbuf_size / (cb_data->Nds * sample_bytesize); - const size_t pol_stride = spectra_stride / cb_data->guppiraw_header.metadata.datashape.n_pol; - const size_t ant_stride = pol_stride / cb_data->Nant; - const size_t chan_stride = pol_stride / cb_data->Nf; // Nf is OBSNCHAN if( - guppiraw_write_block_arbitrary( + guppiraw_write_block( cb_data->fd_ics, &cb_data->guppiraw_header_ics, - cb_data->h_pwrbuf, - ant_stride*sample_bytesize, - chan_stride*sample_bytesize, - spectra_stride*sample_bytesize, - pol_stride*sample_bytesize - ) < (int64_t) cb_data->h_pwrbuf_size / (cb_data->per_ant_out ? 1 : cb_data->Nant) + cb_data->h_icsbuf + ) < (int64_t) cb_data->h_pwrbuf_size / cb_data->Nant ) { cb_data->exit_soon = 1; cb_data->output_thread_valid = 0; diff --git a/rawspec_gpu.cu b/rawspec_gpu.cu index 182a28d..ea19be6 100644 --- a/rawspec_gpu.cu +++ b/rawspec_gpu.cu @@ -533,7 +533,7 @@ int rawspec_initialize(rawspec_context * ctx) ctx->Npolout[i] = 1; } } - ctx->complex_output = min(ctx->complex_output, 2); + ctx->complex_output = min(ctx->complex_output, 1); // Validate Ntpb if(ctx->Ntpb == 0) { @@ -805,6 +805,7 @@ int rawspec_initialize(rawspec_context * ctx) ctx->h_pwrbuf_size[i] = abs(ctx->Npolout[i]) * ctx->Nds[i]*ctx->Nts[i]*ctx->Nc*sizeof(float); ctx->h_pwrbuf_size[i] *= 1 + ctx->complex_output; // complex-output conditional double + printf("FFT Host dump buffer[%d] size == %lu\n", i, ctx->h_pwrbuf_size[i]); #ifdef VERBOSE_ALLOC printf("FFT Host dump buffer[%d] size == %lu\n", i, ctx->h_pwrbuf_size[i]); #endif @@ -1589,6 +1590,7 @@ int rawspec_start_processing(rawspec_context * ctx, int fft_dir) int i; int p; int d; + int b, f; float * dst; size_t dpitch; float * src; @@ -1608,6 +1610,9 @@ int rawspec_start_processing(rawspec_context * ctx, int fft_dir) const char complexity_factor = 1 + ctx->complex_output; bool is_last_channel_batch; + printf("Nss: %d\n", gpu_ctx->Nss[0]); + printf("Nas: %d\n", ctx->Nas[0]); + grid_full_pwr.x = ctx->Nb*ctx->Ntpb*ctx->Nbc; grid_full_pwr.z = complexity_factor; @@ -1616,6 +1621,21 @@ int rawspec_start_processing(rawspec_context * ctx, int fft_dir) // For each output product for(i=0; i < ctx->No; i++) { + if(ctx->complex_output) { + // h_pwrbuf shape: slowest[Channels, Spectra, Polarization]fastest + // have to do many 2D copies to achieve this as fine-channels are no longer contiguous + dpitch = ctx->Nds[i] * abs(ctx->Npolout[i]) * sizeof(float) * complexity_factor; + + spitch = sizeof(float) * complexity_factor; + width = sizeof(float) * complexity_factor; + height = ctx->Nts[i]; + } else { + dpitch = ctx->Nts[i] * sizeof(float) * complexity_factor; + + spitch = gpu_ctx->Nss[i] * ctx->Nts[i] * sizeof(float) * complexity_factor; + height = ctx->Nbc; + } + // Length of an FFT output buffer when abs(Npotout)==4, must be 0 when // Npolout==1 fft_outbuf_length = ctx->Npolout[i] == 1 ? 0 : ctx->Nb*ctx->Ntpb*ctx->Nbc; @@ -1663,7 +1683,7 @@ int rawspec_start_processing(rawspec_context * ctx, int fft_dir) gpu_ctx->nthreads[i], 0, gpu_ctx->compute_stream>>> ( - gpu_ctx->d_pwr_out[i] + p*ctx->Nb*ctx->Ntpb*ctx->Nbc*complexity_factor, + gpu_ctx->d_pwr_out[i] + p*fft_outbuf_length*complexity_factor, MIN(ctx->Nas[i], gpu_ctx->Nss[i]), // Na complexity_factor*ctx->Nts[i], // xpitch complexity_factor*ctx->Nas[i]*ctx->Nts[i], // ypitch @@ -1720,33 +1740,69 @@ int rawspec_start_processing(rawspec_context * ctx, int fft_dir) // Copy integrated power spectra (or spectrum) to host. This is done as // two 2D copies to get channel 0 in the center of the spectrum. Special // care is taken in the unlikely event that Nt is odd. - src = gpu_ctx->d_pwr_out[i] + (p*ctx->Nb*ctx->Ntpb*ctx->Nbc)*complexity_factor; - dst = ctx->h_pwrbuf[i] + ((p*ctx->Nts[i]*ctx->Nc) + (c*ctx->Nts[i]))*complexity_factor; - dpitch = ctx->Nts[i] * sizeof(float) * complexity_factor; - spitch = gpu_ctx->Nss[i] * dpitch; - height = ctx->Nbc; + src = gpu_ctx->d_pwr_out[i] + (p*fft_outbuf_length)*complexity_factor; + if(ctx->complex_output) { + dst = ctx->h_pwrbuf[i] + (c*ctx->Nts[i]*ctx->Nds[i]*abs(ctx->Npolout[i]) + p)*complexity_factor; + } else { + dst = ctx->h_pwrbuf[i] + ((p*ctx->Nts[i]*ctx->Nc) + (c*ctx->Nts[i]))*complexity_factor; + } for(d=0; d < ctx->Nds[i]; d++) { if(ctx->complex_output) { // GUPPI RAW output, don't translate channels - width = dpitch; - cuda_rc = cudaMemcpy2DAsync(dst, - dpitch, - src, - spitch, - width, - height, - cudaMemcpyDeviceToHost, - gpu_ctx->compute_stream); - - if(cuda_rc != cudaSuccess) { - PRINT_CUDA_ERRMSG(cuda_rc); - rawspec_cleanup(ctx); - return 1; + // h_pwrbuf shape: slowest[Channels, Spectra, Polarization]fastest + for(b=0; b < ctx->Nbc; b++) { + #if 0 + for(f=0; f < ctx->Nts[i]; f++) { + cuda_rc = cudaMemcpyAsync( + dst + ((b * ctx->Nts[i] + f) * ctx->Nds[i] * abs(ctx->Npolout[i]) * complexity_factor), + src + ((b*ctx->Nts[i]*gpu_ctx->Nss[i] + f) * complexity_factor), + complexity_factor*sizeof(float), + cudaMemcpyDeviceToHost, + gpu_ctx->compute_stream + ); + // cuda_rc = cudaMemsetAsync( + // dst + ((b * ctx->Nts[i] + f) * ctx->Nds[i] * abs(ctx->Npolout[i]) * complexity_factor), + // 0, + // complexity_factor*sizeof(float), + // gpu_ctx->compute_stream + // ); + // cuda_rc = cudaMemsetAsync( + // dst + ((b * ctx->Nts[i] + f) * ctx->Nds[i] * abs(ctx->Npolout[i]) * complexity_factor), + // b+p, + // 1, + // gpu_ctx->compute_stream + // ); + + if(cuda_rc != cudaSuccess) { + PRINT_CUDA_ERRMSG(cuda_rc); + rawspec_cleanup(ctx); + return 1; + } + } + #else + cuda_rc = cudaMemcpy2DAsync( + dst + (b * ctx->Nts[i] * ctx->Nds[i] * abs(ctx->Npolout[i]) * complexity_factor), + dpitch, + src + (b*gpu_ctx->Nss[i]*ctx->Nts[i]*complexity_factor), + spitch, + width, + height, + cudaMemcpyDeviceToHost, + gpu_ctx->compute_stream + ); + + if(cuda_rc != cudaSuccess) { + PRINT_CUDA_ERRMSG(cuda_rc); + rawspec_cleanup(ctx); + return 1; + } + #endif } } else { + // h_pwrbuf shape: slowest[Spectra, Polarization, Channels]fastest // Lo to hi width = ((ctx->Nts[i]+1) / 2) * sizeof(float); cuda_rc = cudaMemcpy2DAsync(dst + ctx->Nts[i]/2, @@ -1782,9 +1838,13 @@ int rawspec_start_processing(rawspec_context * ctx, int fft_dir) } } - // Increment src and dst pointers + // Increment src and dst pointers, striding spectra src += ctx->Nts[i] * ctx->Nas[i] * complexity_factor; - dst += abs(ctx->Npolout[i]) * ctx->Nts[i] * ctx->Nc * complexity_factor; + if(ctx->complex_output) { + dst += abs(ctx->Npolout[i]) * complexity_factor; + } else { + dst += abs(ctx->Npolout[i]) * ctx->Nts[i] * ctx->Nc * complexity_factor; + } } } @@ -1802,7 +1862,7 @@ int rawspec_start_processing(rawspec_context * ctx, int fft_dir) // Add power buffer clearing cudaMemset call to stream cuda_rc = cudaMemsetAsync(gpu_ctx->d_pwr_out[i], 0, - abs(ctx->Npolout[i])*ctx->Nb*ctx->Ntpb*ctx->Nbc*sizeof(float), + abs(ctx->Npolout[i])*ctx->Nb*ctx->Ntpb*ctx->Nbc*sizeof(float)*complexity_factor, gpu_ctx->compute_stream); if(cuda_rc != cudaSuccess) { @@ -1812,7 +1872,7 @@ int rawspec_start_processing(rawspec_context * ctx, int fft_dir) if(is_last_channel_batch && ctx->incoherently_sum){ // Add ics buffer clearing cudaMemset call to stream cuda_rc = cudaMemsetAsync(gpu_ctx->d_ics_out[i], 0, - abs(ctx->Npolout[i])*ctx->Nb*ctx->Ntpb*ctx->Nc*sizeof(float)/ctx->Nant, + abs(ctx->Npolout[i])*ctx->Nb*ctx->Ntpb*ctx->Nc*sizeof(float)*complexity_factor/ctx->Nant, gpu_ctx->compute_stream); if(cuda_rc != cudaSuccess) { @@ -1826,7 +1886,7 @@ int rawspec_start_processing(rawspec_context * ctx, int fft_dir) // Cache d_pwr_buf to d_prev_pwr_out_cache for later // mem2dAccumMoveFloat memset(0)s the source (d_pwr_buf) mem2dAccumMoveFloat<<compute_stream - >>>(gpu_ctx->d_prev_pwr_out_cache[i] + (c * abs(ctx->Npolout[i]) * ctx->Nb*ctx->Ntpb), + >>>(gpu_ctx->d_prev_pwr_out_cache[i] + (c * abs(ctx->Npolout[i]) * ctx->Nb*ctx->Ntpb)*complexity_factor, ctx->Nb*ctx->Ntpb*ctx->Nc, gpu_ctx->d_pwr_out[i], ctx->Nb*ctx->Ntpb*ctx->Nbc diff --git a/rawspec_rawutils.c b/rawspec_rawutils.c index 3b3f332..ac7c47b 100644 --- a/rawspec_rawutils.c +++ b/rawspec_rawutils.c @@ -173,7 +173,7 @@ void _rawspec_header_parse_metadata(const char* entry, void* raw_hdr_void) { rawspec_raw_hdr_t * raw_hdr = (rawspec_raw_hdr_t*) raw_hdr_void; if(((uint64_t*)entry)[0] == KEY_UINT64_BLOCSIZE) - hgeti4(entry, "BLOCSIZE", &raw_hdr->blocsize); + hgetu8(entry, "BLOCSIZE", &raw_hdr->blocsize); else if(((uint64_t*)entry)[0] == KEY_UINT64_NPOL) hgeti4(entry, "NPOL", &raw_hdr->npol); else if(((uint64_t*)entry)[0] == KEY_UINT64_OBSNCHAN)