zzh
/
fast-bwa
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Reduce dependency on utils.h - new malloc wrapping scheme. 

Remove xmalloc, xcalloc, xrealloc and xstrdup from utils.h and revert calls
to the normal malloc, calloc, realloc, strdup.  Add new files malloc_wrap.[ch]
with the wrapper functions.  malloc_wrap.h #defines malloc etc. to the
wrapper, but only if USE_MALLOC_WRAPPERS has been defined.

Put #include "malloc_wrap.h" in any file that uses *alloc or strdup.  This
is also in a #ifdef USE_MALLOC_WRAPPERS ... #endif block to make using the
wrappers optional.  Add -DUSE_MALLOC_WRAPPERS into the makefile so they
should normally get added.

This is an improvement on the previous method as we now don't need to
worry about stray function calls that were not changed to the wrapped version
and the code will still work even if the wrapping is disabled.

Other possible methods of doing this are using malloc_hook (glibc-specific),
adding -include malloc_wrap.h to the gcc command-line (somewhat
gcc-specific) or making our own malloc function and using dlopen (scary).
This way is probably the most portable.
This commit is contained in:
Rob Davies 2013-05-02 15:12:01 +01:00
parent 0aa7e0a402
commit 96e445d9e4
37 changed files with 516 additions and 342 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

70
Makefile
View File

@ -1,12 +1,13 @@
CC= gcc CC= gcc
CFLAGS= -g -Wall -O2 CFLAGS= -g -Wall -O2
WRAP_MALLOC= -DUSE_MALLOC_WRAPPERS
AR= ar AR= ar
DFLAGS= -DHAVE_PTHREAD DFLAGS= -DHAVE_PTHREAD $(WRAP_MALLOC)
LOBJS= utils.o kstring.o ksw.o bwt.o bntseq.o bwa.o bwamem.o bwamem_pair.o LOBJS= utils.o kstring.o ksw.o bwt.o bntseq.o bwa.o bwamem.o bwamem_pair.o
AOBJS= QSufSort.o bwt_gen.o bwase.o bwaseqio.o bwtgap.o bwtaln.o bamlite.o \ AOBJS= QSufSort.o bwt_gen.o bwase.o bwaseqio.o bwtgap.o bwtaln.o bamlite.o \
is.o bwtindex.o bwape.o kopen.o pemerge.o \ is.o bwtindex.o bwape.o kopen.o pemerge.o \
bwtsw2_core.o bwtsw2_main.o bwtsw2_aux.o bwt_lite.o \ bwtsw2_core.o bwtsw2_main.o bwtsw2_aux.o bwt_lite.o \
bwtsw2_chain.o fastmap.o bwtsw2_pair.o bwtsw2_chain.o fastmap.o bwtsw2_pair.o malloc_wrap.o
PROG= bwa PROG= bwa
INCLUDES= INCLUDES=
LIBS= -lm -lz -lpthread LIBS= -lm -lz -lpthread
@ -32,39 +33,44 @@ clean:
rm -f gmon.out *.o a.out $(PROG) *~ *.a rm -f gmon.out *.o a.out $(PROG) *~ *.a
depend: depend:
( LC_ALL=C ; export LC_ALL; makedepend -Y -- $(CFLAGS) -- *.c ) ( LC_ALL=C ; export LC_ALL; makedepend -Y -- $(CFLAGS) $(DFLAGS) -- *.c )
# DO NOT DELETE THIS LINE -- make depend depends on it. # DO NOT DELETE THIS LINE -- make depend depends on it.
QSufSort.o: QSufSort.h QSufSort.o: QSufSort.h
bamlite.o: utils.h bamlite.h bamlite.o: bamlite.h utils.h malloc_wrap.h
bntseq.o: bntseq.h utils.h kseq.h bntseq.o: bntseq.h utils.h kseq.h malloc_wrap.h
bwa.o: bntseq.h bwa.h bwt.h ksw.h utils.h kseq.h bwa.o: bntseq.h bwa.h bwt.h ksw.h malloc_wrap.h utils.h kseq.h
bwamem.o: kstring.h utils.h bwamem.h bwt.h bntseq.h bwa.h ksw.h kvec.h bwamem.o: kstring.h malloc_wrap.h bwamem.h bwt.h bntseq.h bwa.h ksw.h kvec.h
bwamem.o: ksort.h kbtree.h bwamem.o: ksort.h utils.h kbtree.h
bwamem_pair.o: kstring.h utils.h bwamem.h bwt.h bntseq.h bwa.h kvec.h ksw.h bwamem_pair.o: kstring.h malloc_wrap.h bwamem.h bwt.h bntseq.h bwa.h kvec.h
bwape.o: bwtaln.h bwt.h kvec.h bntseq.h utils.h bwase.h bwa.h ksw.h khash.h bwamem_pair.o: utils.h ksw.h
bwase.o: bwase.h bntseq.h bwt.h bwtaln.h utils.h kstring.h bwa.h ksw.h bwape.o: bwtaln.h bwt.h kvec.h malloc_wrap.h bntseq.h utils.h bwase.h bwa.h
bwaseqio.o: bwtaln.h bwt.h utils.h bamlite.h kseq.h bwape.o: ksw.h khash.h
bwt.o: utils.h bwt.h kvec.h bwase.o: bwase.h bntseq.h bwt.h bwtaln.h utils.h kstring.h malloc_wrap.h
bwt_gen.o: QSufSort.h utils.h bwase.o: bwa.h ksw.h
bwt_lite.o: bwt_lite.h utils.h bwaseqio.o: bwtaln.h bwt.h utils.h bamlite.h malloc_wrap.h kseq.h
bwtaln.o: bwtaln.h bwt.h bwtgap.h utils.h bwa.h bntseq.h bwt.o: utils.h bwt.h kvec.h malloc_wrap.h
bwtgap.o: bwtgap.h bwt.h bwtaln.h utils.h bwt_gen.o: QSufSort.h utils.h malloc_wrap.h
bwtindex.o: bntseq.h bwt.h utils.h bwt_lite.o: bwt_lite.h malloc_wrap.h
bwtsw2_aux.o: bntseq.h bwt_lite.h utils.h bwtsw2.h bwt.h kstring.h bwa.h bwtaln.o: bwtaln.h bwt.h bwtgap.h utils.h bwa.h bntseq.h malloc_wrap.h
bwtsw2_aux.o: ksw.h kseq.h ksort.h bwtgap.o: bwtgap.h bwt.h bwtaln.h malloc_wrap.h
bwtsw2_chain.o: bwtsw2.h bntseq.h bwt_lite.h bwt.h utils.h ksort.h bwtindex.o: bntseq.h bwt.h utils.h malloc_wrap.h
bwtsw2_core.o: bwt_lite.h bwtsw2.h bntseq.h bwt.h kvec.h utils.h khash.h bwtsw2_aux.o: bntseq.h bwt_lite.h utils.h bwtsw2.h bwt.h kstring.h
bwtsw2_core.o: ksort.h bwtsw2_aux.o: malloc_wrap.h bwa.h ksw.h kseq.h ksort.h
bwtsw2_chain.o: bwtsw2.h bntseq.h bwt_lite.h bwt.h malloc_wrap.h ksort.h
bwtsw2_core.o: bwt_lite.h bwtsw2.h bntseq.h bwt.h kvec.h malloc_wrap.h
bwtsw2_core.o: khash.h ksort.h
bwtsw2_main.o: bwt.h bwtsw2.h bntseq.h bwt_lite.h utils.h bwa.h bwtsw2_main.o: bwt.h bwtsw2.h bntseq.h bwt_lite.h utils.h bwa.h
bwtsw2_pair.o: utils.h bwt.h bntseq.h bwtsw2.h bwt_lite.h kstring.h ksw.h bwtsw2_pair.o: utils.h bwt.h bntseq.h bwtsw2.h bwt_lite.h kstring.h
example.o: bwamem.h bwt.h bntseq.h bwa.h kseq.h utils.h bwtsw2_pair.o: malloc_wrap.h ksw.h
fastmap.o: bwa.h bntseq.h bwt.h bwamem.h kvec.h utils.h kseq.h example.o: bwamem.h bwt.h bntseq.h bwa.h kseq.h malloc_wrap.h
is.o: utils.h fastmap.o: bwa.h bntseq.h bwt.h bwamem.h kvec.h malloc_wrap.h utils.h kseq.h
kopen.o: utils.h is.o: utils.h malloc_wrap.h
kstring.o: kstring.h utils.h kopen.o: malloc_wrap.h
ksw.o: ksw.h utils.h kstring.o: kstring.h malloc_wrap.h
ksw.o: ksw.h malloc_wrap.h
main.o: utils.h main.o: utils.h
pemerge.o: ksw.h kseq.h utils.h kstring.h bwa.h bntseq.h bwt.h malloc_wrap.o: malloc_wrap.h
utils.o: utils.h ksort.h kseq.h pemerge.o: ksw.h malloc_wrap.h kseq.h kstring.h bwa.h bntseq.h bwt.h utils.h
utils.o: utils.h ksort.h malloc_wrap.h kseq.h

17
bamlite.c
View File

@ -2,9 +2,12 @@
#include <ctype.h> #include <ctype.h>
#include <string.h> #include <string.h>
#include <stdio.h> #include <stdio.h>
#include "utils.h"
#include "bamlite.h" #include "bamlite.h"
#ifdef USE_MALLOC_WRAPPERS
# include "malloc_wrap.h"
#endif
/********************* /*********************
* from bam_endian.c * * from bam_endian.c *
*********************/ *********************/
@ -54,7 +57,7 @@ int bam_is_be;
bam_header_t *bam_header_init() bam_header_t *bam_header_init()
{ {
bam_is_be = bam_is_big_endian(); bam_is_be = bam_is_big_endian();
return (bam_header_t*)xcalloc(1, sizeof(bam_header_t)); return (bam_header_t*)calloc(1, sizeof(bam_header_t));
} }
void bam_header_destroy(bam_header_t *header) void bam_header_destroy(bam_header_t *header)
@ -87,17 +90,17 @@ bam_header_t *bam_header_read(bamFile fp)
// read plain text and the number of reference sequences // read plain text and the number of reference sequences
if (bam_read(fp, &header->l_text, 4) != 4) goto fail; if (bam_read(fp, &header->l_text, 4) != 4) goto fail;
if (bam_is_be) bam_swap_endian_4p(&header->l_text); if (bam_is_be) bam_swap_endian_4p(&header->l_text);
header->text = (char*)xcalloc(header->l_text + 1, 1); header->text = (char*)calloc(header->l_text + 1, 1);
if (bam_read(fp, header->text, header->l_text) != header->l_text) goto fail; if (bam_read(fp, header->text, header->l_text) != header->l_text) goto fail;
if (bam_read(fp, &header->n_targets, 4) != 4) goto fail; if (bam_read(fp, &header->n_targets, 4) != 4) goto fail;
if (bam_is_be) bam_swap_endian_4p(&header->n_targets); if (bam_is_be) bam_swap_endian_4p(&header->n_targets);
// read reference sequence names and lengths // read reference sequence names and lengths
header->target_name = (char**)xcalloc(header->n_targets, sizeof(char*)); header->target_name = (char**)calloc(header->n_targets, sizeof(char*));
header->target_len = (uint32_t*)xcalloc(header->n_targets, 4); header->target_len = (uint32_t*)calloc(header->n_targets, 4);
for (i = 0; i != header->n_targets; ++i) { for (i = 0; i != header->n_targets; ++i) {
if (bam_read(fp, &name_len, 4) != 4) goto fail; if (bam_read(fp, &name_len, 4) != 4) goto fail;
if (bam_is_be) bam_swap_endian_4p(&name_len); if (bam_is_be) bam_swap_endian_4p(&name_len);
header->target_name[i] = (char*)xcalloc(name_len, 1); header->target_name[i] = (char*)calloc(name_len, 1);
if (bam_read(fp, header->target_name[i], name_len) != name_len) { if (bam_read(fp, header->target_name[i], name_len) != name_len) {
goto fail; goto fail;
} }
@ -152,7 +155,7 @@ int bam_read1(bamFile fp, bam1_t *b)
if (b->m_data < b->data_len) { if (b->m_data < b->data_len) {
b->m_data = b->data_len; b->m_data = b->data_len;
kroundup32(b->m_data); kroundup32(b->m_data);
b->data = (uint8_t*)xrealloc(b->data, b->m_data); b->data = (uint8_t*)realloc(b->data, b->m_data);
} }
if (bam_read(fp, b->data, b->data_len) != b->data_len) return -4; if (bam_read(fp, b->data, b->data_len) != b->data_len) return -4;
b->l_aux = b->data_len - c->n_cigar * 4 - c->l_qname - c->l_qseq - (c->l_qseq+1)/2; b->l_aux = b->data_len - c->n_cigar * 4 - c->l_qname - c->l_qseq - (c->l_qseq+1)/2;

6
bamlite.h
View File

@ -5,6 +5,10 @@
#include <zlib.h> #include <zlib.h>
#include "utils.h" #include "utils.h"
#ifdef USE_MALLOC_WRAPPERS
# include "malloc_wrap.h"
#endif
typedef gzFile bamFile; typedef gzFile bamFile;
#define bam_open(fn, mode) xzopen(fn, mode) #define bam_open(fn, mode) xzopen(fn, mode)
#define bam_dopen(fd, mode) gzdopen(fd, mode) #define bam_dopen(fd, mode) gzdopen(fd, mode)
@ -72,7 +76,7 @@ typedef struct {
#define bam1_seqi(s, i) ((s)[(i)/2] >> 4*(1-(i)%2) & 0xf) #define bam1_seqi(s, i) ((s)[(i)/2] >> 4*(1-(i)%2) & 0xf)
#define bam1_aux(b) ((b)->data + (b)->core.n_cigar*4 + (b)->core.l_qname + (b)->core.l_qseq + ((b)->core.l_qseq + 1)/2) #define bam1_aux(b) ((b)->data + (b)->core.n_cigar*4 + (b)->core.l_qname + (b)->core.l_qseq + ((b)->core.l_qseq + 1)/2)
#define bam_init1() ((bam1_t*)xcalloc(1, sizeof(bam1_t))) #define bam_init1() ((bam1_t*)calloc(1, sizeof(bam1_t)))
#define bam_destroy1(b) do { \ #define bam_destroy1(b) do { \
if (b) { free((b)->data); free(b); } \ if (b) { free((b)->data); free(b); } \
} while (0) } while (0)

38
bntseq.c
View File

@ -37,6 +37,10 @@
#include "kseq.h" #include "kseq.h"
KSEQ_DECLARE(gzFile) KSEQ_DECLARE(gzFile)
#ifdef USE_MALLOC_WRAPPERS
# include "malloc_wrap.h"
#endif
unsigned char nst_nt4_table[256] = { unsigned char nst_nt4_table[256] = {
4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,
4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,
@ -97,13 +101,13 @@ bntseq_t *bns_restore_core(const char *ann_filename, const char* amb_filename, c
long long xx; long long xx;
int i; int i;
int scanres; int scanres;
bns = (bntseq_t*)xcalloc(1, sizeof(bntseq_t)); bns = (bntseq_t*)calloc(1, sizeof(bntseq_t));
{ // read .ann { // read .ann
fp = xopen(fname = ann_filename, "r"); fp = xopen(fname = ann_filename, "r");
scanres = fscanf(fp, "%lld%d%u", &xx, &bns->n_seqs, &bns->seed); scanres = fscanf(fp, "%lld%d%u", &xx, &bns->n_seqs, &bns->seed);
if (scanres != 3) goto badread; if (scanres != 3) goto badread;
bns->l_pac = xx; bns->l_pac = xx;
bns->anns = (bntann1_t*)xcalloc(bns->n_seqs, sizeof(bntann1_t)); bns->anns = (bntann1_t*)calloc(bns->n_seqs, sizeof(bntann1_t));
for (i = 0; i < bns->n_seqs; ++i) { for (i = 0; i < bns->n_seqs; ++i) {
bntann1_t *p = bns->anns + i; bntann1_t *p = bns->anns + i;
char *q = str; char *q = str;
@ -111,7 +115,7 @@ bntseq_t *bns_restore_core(const char *ann_filename, const char* amb_filename, c
// read gi and sequence name // read gi and sequence name
scanres = fscanf(fp, "%u%s", &p->gi, str); scanres = fscanf(fp, "%u%s", &p->gi, str);
if (scanres != 2) goto badread; if (scanres != 2) goto badread;
p->name = xstrdup(str); p->name = strdup(str);
// read fasta comments // read fasta comments
while (str - q < sizeof(str) - 1 && (c = fgetc(fp)) != '\n' && c != EOF) *q++ = c; while (str - q < sizeof(str) - 1 && (c = fgetc(fp)) != '\n' && c != EOF) *q++ = c;
while (c != '\n' && c != EOF) c = fgetc(fp); while (c != '\n' && c != EOF) c = fgetc(fp);
@ -120,8 +124,8 @@ bntseq_t *bns_restore_core(const char *ann_filename, const char* amb_filename, c
goto badread; goto badread;
} }
*q = 0; *q = 0;
if (q - str > 1) p->anno = xstrdup(str + 1); // skip leading space if (q - str > 1) p->anno = strdup(str + 1); // skip leading space
else p->anno = xstrdup(""); else p->anno = strdup("");
// read the rest // read the rest
scanres = fscanf(fp, "%lld%d%d", &xx, &p->len, &p->n_ambs); scanres = fscanf(fp, "%lld%d%d", &xx, &p->len, &p->n_ambs);
if (scanres != 3) goto badread; if (scanres != 3) goto badread;
@ -137,7 +141,7 @@ bntseq_t *bns_restore_core(const char *ann_filename, const char* amb_filename, c
if (scanres != 3) goto badread; if (scanres != 3) goto badread;
l_pac = xx; l_pac = xx;
xassert(l_pac == bns->l_pac && n_seqs == bns->n_seqs, "inconsistent .ann and .amb files."); xassert(l_pac == bns->l_pac && n_seqs == bns->n_seqs, "inconsistent .ann and .amb files.");
bns->ambs = bns->n_holes? (bntamb1_t*)xcalloc(bns->n_holes, sizeof(bntamb1_t)) : 0; bns->ambs = bns->n_holes? (bntamb1_t*)calloc(bns->n_holes, sizeof(bntamb1_t)) : 0;
for (i = 0; i < bns->n_holes; ++i) { for (i = 0; i < bns->n_holes; ++i) {
bntamb1_t *p = bns->ambs + i; bntamb1_t *p = bns->ambs + i;
scanres = fscanf(fp, "%lld%d%s", &xx, &p->len, str); scanres = fscanf(fp, "%lld%d%s", &xx, &p->len, str);
@ -193,11 +197,11 @@ static uint8_t *add1(const kseq_t *seq, bntseq_t *bns, uint8_t *pac, int64_t *m_
int i, lasts; int i, lasts;
if (bns->n_seqs == *m_seqs) { if (bns->n_seqs == *m_seqs) {
*m_seqs <<= 1; *m_seqs <<= 1;
bns->anns = (bntann1_t*)xrealloc(bns->anns, *m_seqs * sizeof(bntann1_t)); bns->anns = (bntann1_t*)realloc(bns->anns, *m_seqs * sizeof(bntann1_t));
} }
p = bns->anns + bns->n_seqs; p = bns->anns + bns->n_seqs;
p->name = xstrdup((char*)seq->name.s); p->name = strdup((char*)seq->name.s);
p->anno = seq->comment.s? xstrdup((char*)seq->comment.s) : xstrdup("(null)"); p->anno = seq->comment.s? strdup((char*)seq->comment.s) : strdup("(null)");
p->gi = 0; p->len = seq->seq.l; p->gi = 0; p->len = seq->seq.l;
p->offset = (bns->n_seqs == 0)? 0 : (p-1)->offset + (p-1)->len; p->offset = (bns->n_seqs == 0)? 0 : (p-1)->offset + (p-1)->len;
p->n_ambs = 0; p->n_ambs = 0;
@ -209,7 +213,7 @@ static uint8_t *add1(const kseq_t *seq, bntseq_t *bns, uint8_t *pac, int64_t *m_
} else { } else {
if (bns->n_holes == *m_holes) { if (bns->n_holes == *m_holes) {
(*m_holes) <<= 1; (*m_holes) <<= 1;
bns->ambs = (bntamb1_t*)xrealloc(bns->ambs, (*m_holes) * sizeof(bntamb1_t)); bns->ambs = (bntamb1_t*)realloc(bns->ambs, (*m_holes) * sizeof(bntamb1_t));
} }
*q = bns->ambs + bns->n_holes; *q = bns->ambs + bns->n_holes;
(*q)->len = 1; (*q)->len = 1;
@ -224,7 +228,7 @@ static uint8_t *add1(const kseq_t *seq, bntseq_t *bns, uint8_t *pac, int64_t *m_
if (c >= 4) c = lrand48()&3; if (c >= 4) c = lrand48()&3;
if (bns->l_pac == *m_pac) { // double the pac size if (bns->l_pac == *m_pac) { // double the pac size
*m_pac <<= 1; *m_pac <<= 1;
pac = xrealloc(pac, *m_pac/4); pac = realloc(pac, *m_pac/4);
memset(pac + bns->l_pac/4, 0, (*m_pac - bns->l_pac)/4); memset(pac + bns->l_pac/4, 0, (*m_pac - bns->l_pac)/4);
} }
_set_pac(pac, bns->l_pac, c); _set_pac(pac, bns->l_pac, c);
@ -249,13 +253,13 @@ int64_t bns_fasta2bntseq(gzFile fp_fa, const char *prefix, int for_only)
// initialization // initialization
seq = kseq_init(fp_fa); seq = kseq_init(fp_fa);
bns = (bntseq_t*)xcalloc(1, sizeof(bntseq_t)); bns = (bntseq_t*)calloc(1, sizeof(bntseq_t));
bns->seed = 11; // fixed seed for random generator bns->seed = 11; // fixed seed for random generator
srand48(bns->seed); srand48(bns->seed);
m_seqs = m_holes = 8; m_pac = 0x10000; m_seqs = m_holes = 8; m_pac = 0x10000;
bns->anns = (bntann1_t*)xcalloc(m_seqs, sizeof(bntann1_t)); bns->anns = (bntann1_t*)calloc(m_seqs, sizeof(bntann1_t));
bns->ambs = (bntamb1_t*)xcalloc(m_holes, sizeof(bntamb1_t)); bns->ambs = (bntamb1_t*)calloc(m_holes, sizeof(bntamb1_t));
pac = xcalloc(m_pac/4, 1); pac = calloc(m_pac/4, 1);
q = bns->ambs; q = bns->ambs;
strcpy(name, prefix); strcat(name, ".pac"); strcpy(name, prefix); strcat(name, ".pac");
fp = xopen(name, "wb"); fp = xopen(name, "wb");
@ -263,7 +267,7 @@ int64_t bns_fasta2bntseq(gzFile fp_fa, const char *prefix, int for_only)
while (kseq_read(seq) >= 0) pac = add1(seq, bns, pac, &m_pac, &m_seqs, &m_holes, &q); while (kseq_read(seq) >= 0) pac = add1(seq, bns, pac, &m_pac, &m_seqs, &m_holes, &q);
if (!for_only) { // add the reverse complemented sequence if (!for_only) { // add the reverse complemented sequence
m_pac = (bns->l_pac * 2 + 3) / 4 * 4; m_pac = (bns->l_pac * 2 + 3) / 4 * 4;
pac = xrealloc(pac, m_pac/4); pac = realloc(pac, m_pac/4);
memset(pac + (bns->l_pac+3)/4, 0, (m_pac - (bns->l_pac+3)/4*4) / 4); memset(pac + (bns->l_pac+3)/4, 0, (m_pac - (bns->l_pac+3)/4*4) / 4);
for (l = bns->l_pac - 1; l >= 0; --l, ++bns->l_pac) for (l = bns->l_pac - 1; l >= 0; --l, ++bns->l_pac)
_set_pac(pac, bns->l_pac, 3-_get_pac(pac, l)); _set_pac(pac, bns->l_pac, 3-_get_pac(pac, l));
@ -357,7 +361,7 @@ uint8_t *bns_get_seq(int64_t l_pac, const uint8_t *pac, int64_t beg, int64_t end
if (beg >= l_pac || end <= l_pac) { if (beg >= l_pac || end <= l_pac) {
int64_t k, l = 0; int64_t k, l = 0;
*len = end - beg; *len = end - beg;
seq = xmalloc(end - beg); seq = malloc(end - beg);
if (beg >= l_pac) { // reverse strand if (beg >= l_pac) { // reverse strand
int64_t beg_f = (l_pac<<1) - 1 - end; int64_t beg_f = (l_pac<<1) - 1 - end;
int64_t end_f = (l_pac<<1) - 1 - beg; int64_t end_f = (l_pac<<1) - 1 - beg;

26
bwa.c
View File

@ -7,6 +7,10 @@
#include "ksw.h" #include "ksw.h"
#include "utils.h" #include "utils.h"
#ifdef USE_MALLOC_WRAPPERS
# include "malloc_wrap.h"
#endif
int bwa_verbose = 3; int bwa_verbose = 3;
char bwa_rg_id[256]; char bwa_rg_id[256];
@ -25,10 +29,10 @@ static inline void trim_readno(kstring_t *s)
static inline void kseq2bseq1(const kseq_t *ks, bseq1_t *s) static inline void kseq2bseq1(const kseq_t *ks, bseq1_t *s)
{ // TODO: it would be better to allocate one chunk of memory, but probably it does not matter in practice { // TODO: it would be better to allocate one chunk of memory, but probably it does not matter in practice
s->name = xstrdup(ks->name.s); s->name = strdup(ks->name.s);
s->comment = ks->comment.l? xstrdup(ks->comment.s) : 0; s->comment = ks->comment.l? strdup(ks->comment.s) : 0;
s->seq = xstrdup(ks->seq.s); s->seq = strdup(ks->seq.s);
s->qual = ks->qual.l? xstrdup(ks->qual.s) : 0; s->qual = ks->qual.l? strdup(ks->qual.s) : 0;
s->l_seq = strlen(s->seq); s->l_seq = strlen(s->seq);
} }
@ -45,7 +49,7 @@ bseq1_t *bseq_read(int chunk_size, int *n_, void *ks1_, void *ks2_)
} }
if (n >= m) { if (n >= m) {
m = m? m<<1 : 256; m = m? m<<1 : 256;
seqs = xrealloc(seqs, m * sizeof(bseq1_t)); seqs = realloc(seqs, m * sizeof(bseq1_t));
} }
trim_readno(&ks->name); trim_readno(&ks->name);
kseq2bseq1(ks, &seqs[n]); kseq2bseq1(ks, &seqs[n]);
@ -98,7 +102,7 @@ uint32_t *bwa_gen_cigar(const int8_t mat[25], int q, int r, int w_, int64_t l_pa
tmp = rseq[i], rseq[i] = rseq[rlen - 1 - i], rseq[rlen - 1 - i] = tmp; tmp = rseq[i], rseq[i] = rseq[rlen - 1 - i], rseq[rlen - 1 - i] = tmp;
} }
if (l_query == re - rb && w_ == 0) { // no gap; no need to do DP if (l_query == re - rb && w_ == 0) { // no gap; no need to do DP
cigar = xmalloc(4); cigar = malloc(4);
cigar[0] = l_query<<4 | 0; cigar[0] = l_query<<4 | 0;
*n_cigar = 1; *n_cigar = 1;
for (i = 0, *score = 0; i < l_query; ++i) for (i = 0, *score = 0; i < l_query; ++i)
@ -205,7 +209,7 @@ char *bwa_idx_infer_prefix(const char *hint)
int l_hint; int l_hint;
FILE *fp; FILE *fp;
l_hint = strlen(hint); l_hint = strlen(hint);
prefix = xmalloc(l_hint + 3 + 4 + 1); prefix = malloc(l_hint + 3 + 4 + 1);
strcpy(prefix, hint); strcpy(prefix, hint);
strcpy(prefix + l_hint, ".64.bwt"); strcpy(prefix + l_hint, ".64.bwt");
if ((fp = fopen(prefix, "rb")) != 0) { if ((fp = fopen(prefix, "rb")) != 0) {
@ -234,7 +238,7 @@ bwt_t *bwa_idx_load_bwt(const char *hint)
if (bwa_verbose >= 1) fprintf(stderr, "[E::%s] fail to locate the index files\n", __func__); if (bwa_verbose >= 1) fprintf(stderr, "[E::%s] fail to locate the index files\n", __func__);
return 0; return 0;
} }
tmp = xcalloc(strlen(prefix) + 5, 1); tmp = calloc(strlen(prefix) + 5, 1);
strcat(strcpy(tmp, prefix), ".bwt"); // FM-index strcat(strcpy(tmp, prefix), ".bwt"); // FM-index
bwt = bwt_restore_bwt(tmp); bwt = bwt_restore_bwt(tmp);
strcat(strcpy(tmp, prefix), ".sa"); // partial suffix array (SA) strcat(strcpy(tmp, prefix), ".sa"); // partial suffix array (SA)
@ -252,12 +256,12 @@ bwaidx_t *bwa_idx_load(const char *hint, int which)
if (bwa_verbose >= 1) fprintf(stderr, "[E::%s] fail to locate the index files\n", __func__); if (bwa_verbose >= 1) fprintf(stderr, "[E::%s] fail to locate the index files\n", __func__);
return 0; return 0;
} }
idx = xcalloc(1, sizeof(bwaidx_t)); idx = calloc(1, sizeof(bwaidx_t));
if (which & BWA_IDX_BWT) idx->bwt = bwa_idx_load_bwt(hint); if (which & BWA_IDX_BWT) idx->bwt = bwa_idx_load_bwt(hint);
if (which & BWA_IDX_BNS) { if (which & BWA_IDX_BNS) {
idx->bns = bns_restore(prefix); idx->bns = bns_restore(prefix);
if (which & BWA_IDX_PAC) { if (which & BWA_IDX_PAC) {
idx->pac = xcalloc(idx->bns->l_pac/4+1, 1); idx->pac = calloc(idx->bns->l_pac/4+1, 1);
err_fread_noeof(idx->pac, 1, idx->bns->l_pac/4+1, idx->bns->fp_pac); // concatenated 2-bit encoded sequence err_fread_noeof(idx->pac, 1, idx->bns->l_pac/4+1, idx->bns->fp_pac); // concatenated 2-bit encoded sequence
err_fclose(idx->bns->fp_pac); err_fclose(idx->bns->fp_pac);
idx->bns->fp_pac = 0; idx->bns->fp_pac = 0;
@ -312,7 +316,7 @@ char *bwa_set_rg(const char *s)
if (bwa_verbose >= 1) fprintf(stderr, "[E::%s] the read group line is not started with @RG\n", __func__); if (bwa_verbose >= 1) fprintf(stderr, "[E::%s] the read group line is not started with @RG\n", __func__);
goto err_set_rg; goto err_set_rg;
} }
rg_line = xstrdup(s); rg_line = strdup(s);
bwa_escape(rg_line); bwa_escape(rg_line);
if ((p = strstr(rg_line, "\tID:")) == 0) { if ((p = strstr(rg_line, "\tID:")) == 0) {
if (bwa_verbose >= 1) fprintf(stderr, "[E::%s] no ID at the read group line\n", __func__); if (bwa_verbose >= 1) fprintf(stderr, "[E::%s] no ID at the read group line\n", __func__);

42
bwamem.c
View File

@ -15,6 +15,10 @@
#include "ksort.h" #include "ksort.h"
#include "utils.h" #include "utils.h"
#ifdef USE_MALLOC_WRAPPERS
# include "malloc_wrap.h"
#endif
/* Theory on probability and scoring *ungapped* alignment /* Theory on probability and scoring *ungapped* alignment
* *
* s'(a,b) = log[P(b|a)/P(b)] = log[4P(b|a)], assuming uniform base distribution * s'(a,b) = log[P(b|a)/P(b)] = log[4P(b|a)], assuming uniform base distribution
@ -41,7 +45,7 @@
mem_opt_t *mem_opt_init() mem_opt_t *mem_opt_init()
{ {
mem_opt_t *o; mem_opt_t *o;
o = xcalloc(1, sizeof(mem_opt_t)); o = calloc(1, sizeof(mem_opt_t));
o->flag = 0; o->flag = 0;
o->a = 1; o->b = 4; o->q = 6; o->r = 1; o->w = 100; o->a = 1; o->b = 4; o->q = 6; o->r = 1; o->w = 100;
o->T = 30; o->T = 30;
@ -79,12 +83,12 @@ struct __smem_i {
smem_i *smem_itr_init(const bwt_t *bwt) smem_i *smem_itr_init(const bwt_t *bwt)
{ {
smem_i *itr; smem_i *itr;
itr = xcalloc(1, sizeof(smem_i)); itr = calloc(1, sizeof(smem_i));
itr->bwt = bwt; itr->bwt = bwt;
itr->tmpvec[0] = xcalloc(1, sizeof(bwtintv_v)); itr->tmpvec[0] = calloc(1, sizeof(bwtintv_v));
itr->tmpvec[1] = xcalloc(1, sizeof(bwtintv_v)); itr->tmpvec[1] = calloc(1, sizeof(bwtintv_v));
itr->matches = xcalloc(1, sizeof(bwtintv_v)); itr->matches = calloc(1, sizeof(bwtintv_v));
itr->sub = xcalloc(1, sizeof(bwtintv_v)); itr->sub = calloc(1, sizeof(bwtintv_v));
return itr; return itr;
} }
@ -181,7 +185,7 @@ static int test_and_merge(const mem_opt_t *opt, int64_t l_pac, mem_chain_t *c, c
if (y >= 0 && x - y <= opt->w && y - x <= opt->w && x - last->len < opt->max_chain_gap && y - last->len < opt->max_chain_gap) { // grow the chain if (y >= 0 && x - y <= opt->w && y - x <= opt->w && x - last->len < opt->max_chain_gap && y - last->len < opt->max_chain_gap) { // grow the chain
if (c->n == c->m) { if (c->n == c->m) {
c->m <<= 1; c->m <<= 1;
c->seeds = xrealloc(c->seeds, c->m * sizeof(mem_seed_t)); c->seeds = realloc(c->seeds, c->m * sizeof(mem_seed_t));
} }
c->seeds[c->n++] = *p; c->seeds[c->n++] = *p;
return 1; return 1;
@ -215,7 +219,7 @@ static void mem_insert_seed(const mem_opt_t *opt, int64_t l_pac, kbtree_t(chn) *
} else to_add = 1; } else to_add = 1;
if (to_add) { // add the seed as a new chain if (to_add) { // add the seed as a new chain
tmp.n = 1; tmp.m = 4; tmp.n = 1; tmp.m = 4;
tmp.seeds = xcalloc(tmp.m, sizeof(mem_seed_t)); tmp.seeds = calloc(tmp.m, sizeof(mem_seed_t));
tmp.seeds[0] = s; tmp.seeds[0] = s;
kb_putp(chn, tree, &tmp); kb_putp(chn, tree, &tmp);
} }
@ -283,7 +287,7 @@ int mem_chain_flt(const mem_opt_t *opt, int n_chn, mem_chain_t *chains)
flt_aux_t *a; flt_aux_t *a;
int i, j, n; int i, j, n;
if (n_chn <= 1) return n_chn; // no need to filter if (n_chn <= 1) return n_chn; // no need to filter
a = xmalloc(sizeof(flt_aux_t) * n_chn); a = malloc(sizeof(flt_aux_t) * n_chn);
for (i = 0; i < n_chn; ++i) { for (i = 0; i < n_chn; ++i) {
mem_chain_t *c = &chains[i]; mem_chain_t *c = &chains[i];
int64_t end; int64_t end;
@ -309,7 +313,7 @@ int mem_chain_flt(const mem_opt_t *opt, int n_chn, mem_chain_t *chains)
ks_introsort(mem_flt, n_chn, a); ks_introsort(mem_flt, n_chn, a);
{ // reorder chains such that the best chain appears first { // reorder chains such that the best chain appears first
mem_chain_t *swap; mem_chain_t *swap;
swap = xmalloc(sizeof(mem_chain_t) * n_chn); swap = malloc(sizeof(mem_chain_t) * n_chn);
for (i = 0; i < n_chn; ++i) { for (i = 0; i < n_chn; ++i) {
swap[i] = *((mem_chain_t*)a[i].p); swap[i] = *((mem_chain_t*)a[i].p);
a[i].p = &chains[i]; // as we will memcpy() below, a[i].p is changed a[i].p = &chains[i]; // as we will memcpy() below, a[i].p is changed
@ -512,7 +516,7 @@ void mem_chain2aln(const mem_opt_t *opt, int64_t l_pac, const uint8_t *pac, int
rseq = bns_get_seq(l_pac, pac, rmax[0], rmax[1], &rlen); rseq = bns_get_seq(l_pac, pac, rmax[0], rmax[1], &rlen);
assert(rlen == rmax[1] - rmax[0]); assert(rlen == rmax[1] - rmax[0]);
srt = xmalloc(c->n * 8); srt = malloc(c->n * 8);
for (i = 0; i < c->n; ++i) for (i = 0; i < c->n; ++i)
srt[i] = (uint64_t)c->seeds[i].len<<32 | i; srt[i] = (uint64_t)c->seeds[i].len<<32 | i;
ks_introsort_64(c->n, srt); ks_introsort_64(c->n, srt);
@ -560,10 +564,10 @@ void mem_chain2aln(const mem_opt_t *opt, int64_t l_pac, const uint8_t *pac, int
if (s->qbeg) { // left extension if (s->qbeg) { // left extension
uint8_t *rs, *qs; uint8_t *rs, *qs;
int qle, tle, gtle, gscore; int qle, tle, gtle, gscore;
qs = xmalloc(s->qbeg); qs = malloc(s->qbeg);
for (i = 0; i < s->qbeg; ++i) qs[i] = query[s->qbeg - 1 - i]; for (i = 0; i < s->qbeg; ++i) qs[i] = query[s->qbeg - 1 - i];
tmp = s->rbeg - rmax[0]; tmp = s->rbeg - rmax[0];
rs = xmalloc(tmp); rs = malloc(tmp);
for (i = 0; i < tmp; ++i) rs[i] = rseq[tmp - 1 - i]; for (i = 0; i < tmp; ++i) rs[i] = rseq[tmp - 1 - i];
for (i = 0; i < MAX_BAND_TRY; ++i) { for (i = 0; i < MAX_BAND_TRY; ++i) {
int prev = a->score; int prev = a->score;
@ -838,7 +842,7 @@ mem_alnreg_v mem_align1(const mem_opt_t *opt, const bwt_t *bwt, const bntseq_t *
{ // the difference from mem_align1_core() is that this routine: 1) calls mem_mark_primary_se(); 2) does not modify the input sequence { // the difference from mem_align1_core() is that this routine: 1) calls mem_mark_primary_se(); 2) does not modify the input sequence
mem_alnreg_v ar; mem_alnreg_v ar;
char *seq; char *seq;
seq = xmalloc(l_seq); seq = malloc(l_seq);
memcpy(seq, seq_, l_seq); // makes a copy of seq_ memcpy(seq, seq_, l_seq); // makes a copy of seq_
ar = mem_align1_core(opt, bwt, bns, pac, l_seq, seq); ar = mem_align1_core(opt, bwt, bns, pac, l_seq, seq);
mem_mark_primary_se(opt, ar.n, ar.a); mem_mark_primary_se(opt, ar.n, ar.a);
@ -861,7 +865,7 @@ mem_aln_t mem_reg2aln(const mem_opt_t *opt, const bntseq_t *bns, const uint8_t *
} }
qb = ar->qb, qe = ar->qe; qb = ar->qb, qe = ar->qe;
rb = ar->rb, re = ar->re; rb = ar->rb, re = ar->re;
query = xmalloc(l_query); query = malloc(l_query);
for (i = 0; i < l_query; ++i) // convert to the nt4 encoding for (i = 0; i < l_query; ++i) // convert to the nt4 encoding
query[i] = query_[i] < 5? query_[i] : nst_nt4_table[(int)query_[i]]; query[i] = query_[i] < 5? query_[i] : nst_nt4_table[(int)query_[i]];
a.mapq = ar->secondary < 0? mem_approx_mapq_se(opt, ar) : 0; a.mapq = ar->secondary < 0? mem_approx_mapq_se(opt, ar) : 0;
@ -887,7 +891,7 @@ mem_aln_t mem_reg2aln(const mem_opt_t *opt, const bntseq_t *bns, const uint8_t *
int clip5, clip3; int clip5, clip3;
clip5 = is_rev? l_query - qe : qb; clip5 = is_rev? l_query - qe : qb;
clip3 = is_rev? qb : l_query - qe; clip3 = is_rev? qb : l_query - qe;
a.cigar = xrealloc(a.cigar, 4 * (a.n_cigar + 2)); a.cigar = realloc(a.cigar, 4 * (a.n_cigar + 2));
if (clip5) { if (clip5) {
memmove(a.cigar+1, a.cigar, a.n_cigar * 4); memmove(a.cigar+1, a.cigar, a.n_cigar * 4);
a.cigar[0] = clip5<<4 | (opt->flag&MEM_F_HARDCLIP? 4 : 3); a.cigar[0] = clip5<<4 | (opt->flag&MEM_F_HARDCLIP? 4 : 3);
@ -958,8 +962,8 @@ void mem_process_seqs(const mem_opt_t *opt, const bwt_t *bwt, const bntseq_t *bn
mem_alnreg_v *regs; mem_alnreg_v *regs;
mem_pestat_t pes[4]; mem_pestat_t pes[4];
w = xcalloc(opt->n_threads, sizeof(worker_t)); w = calloc(opt->n_threads, sizeof(worker_t));
regs = xmalloc(n * sizeof(mem_alnreg_v)); regs = malloc(n * sizeof(mem_alnreg_v));
for (i = 0; i < opt->n_threads; ++i) { for (i = 0; i < opt->n_threads; ++i) {
worker_t *p = &w[i]; worker_t *p = &w[i];
p->start = i; p->step = opt->n_threads; p->n = n; p->start = i; p->step = opt->n_threads; p->n = n;
@ -980,7 +984,7 @@ void mem_process_seqs(const mem_opt_t *opt, const bwt_t *bwt, const bntseq_t *bn
#ifdef HAVE_PTHREAD #ifdef HAVE_PTHREAD
} else { } else {
pthread_t *tid; pthread_t *tid;
tid = (pthread_t*)xcalloc(opt->n_threads, sizeof(pthread_t)); tid = (pthread_t*)calloc(opt->n_threads, sizeof(pthread_t));
for (i = 0; i < opt->n_threads; ++i) pthread_create(&tid[i], 0, worker1, &w[i]); for (i = 0; i < opt->n_threads; ++i) pthread_create(&tid[i], 0, worker1, &w[i]);
for (i = 0; i < opt->n_threads; ++i) pthread_join(tid[i], 0); for (i = 0; i < opt->n_threads; ++i) pthread_join(tid[i], 0);
if (opt->flag&MEM_F_PE) { if (opt->flag&MEM_F_PE) {

9
bwamem_pair.c
View File

@ -8,6 +8,11 @@
#include "utils.h" #include "utils.h"
#include "ksw.h" #include "ksw.h"
#ifdef USE_MALLOC_WRAPPERS
# include "malloc_wrap.h"
#endif
#define MIN_RATIO 0.8 #define MIN_RATIO 0.8
#define MIN_DIR_CNT 10 #define MIN_DIR_CNT 10
#define MIN_DIR_RATIO 0.05 #define MIN_DIR_RATIO 0.05
@ -121,7 +126,7 @@ int mem_matesw(const mem_opt_t *opt, int64_t l_pac, const uint8_t *pac, const me
is_rev = (r>>1 != (r&1)); // whether to reverse complement the mate is_rev = (r>>1 != (r&1)); // whether to reverse complement the mate
is_larger = !(r>>1); // whether the mate has larger coordinate is_larger = !(r>>1); // whether the mate has larger coordinate
if (is_rev) { if (is_rev) {
rev = xmalloc(l_ms); // this is the reverse complement of $ms rev = malloc(l_ms); // this is the reverse complement of $ms
for (i = 0; i < l_ms; ++i) rev[l_ms - 1 - i] = ms[i] < 4? 3 - ms[i] : 4; for (i = 0; i < l_ms; ++i) rev[l_ms - 1 - i] = ms[i] < 4? 3 - ms[i] : 4;
seq = rev; seq = rev;
} else seq = (uint8_t*)ms; } else seq = (uint8_t*)ms;
@ -294,7 +299,7 @@ int mem_sam_pe(const mem_opt_t *opt, const bntseq_t *bns, const uint8_t *pac, co
// write SAM // write SAM
h[0] = mem_reg2aln(opt, bns, pac, s[0].l_seq, s[0].seq, &a[0].a[z[0]]); h[0].mapq = q_se[0]; h[0].flag |= 0x40 | extra_flag; h[0] = mem_reg2aln(opt, bns, pac, s[0].l_seq, s[0].seq, &a[0].a[z[0]]); h[0].mapq = q_se[0]; h[0].flag |= 0x40 | extra_flag;
h[1] = mem_reg2aln(opt, bns, pac, s[1].l_seq, s[1].seq, &a[1].a[z[1]]); h[1].mapq = q_se[1]; h[1].flag |= 0x80 | extra_flag; h[1] = mem_reg2aln(opt, bns, pac, s[1].l_seq, s[1].seq, &a[1].a[z[1]]); h[1].mapq = q_se[1]; h[1].flag |= 0x80 | extra_flag;
mem_aln2sam(bns, &str, &s[0], 1, &h[0], 0, &h[1]); s[0].sam = xstrdup(str.s); str.l = 0; mem_aln2sam(bns, &str, &s[0], 1, &h[0], 0, &h[1]); s[0].sam = strdup(str.s); str.l = 0;
mem_aln2sam(bns, &str, &s[1], 1, &h[1], 0, &h[0]); s[1].sam = str.s; mem_aln2sam(bns, &str, &s[1], 1, &h[1], 0, &h[0]); s[1].sam = str.s;
free(h[0].cigar); free(h[1].cigar); free(h[0].cigar); free(h[1].cigar);
} else goto no_pairing; } else goto no_pairing;

24
bwape.c
View File

@ -12,6 +12,10 @@
#include "bwa.h" #include "bwa.h"
#include "ksw.h" #include "ksw.h"
#ifdef USE_MALLOC_WRAPPERS
# include "malloc_wrap.h"
#endif
typedef struct { typedef struct {
int n; int n;
bwtint_t *a; bwtint_t *a;
@ -50,7 +54,7 @@ void bwa_print_sam_PG();
pe_opt_t *bwa_init_pe_opt() pe_opt_t *bwa_init_pe_opt()
{ {
pe_opt_t *po; pe_opt_t *po;
po = (pe_opt_t*)xcalloc(1, sizeof(pe_opt_t)); po = (pe_opt_t*)calloc(1, sizeof(pe_opt_t));
po->max_isize = 500; po->max_isize = 500;
po->force_isize = 0; po->force_isize = 0;
po->max_occ = 100000; po->max_occ = 100000;
@ -83,7 +87,7 @@ static int infer_isize(int n_seqs, bwa_seq_t *seqs[2], isize_info_t *ii, double
ii->avg = ii->std = -1.0; ii->avg = ii->std = -1.0;
ii->low = ii->high = ii->high_bayesian = 0; ii->low = ii->high = ii->high_bayesian = 0;
isizes = (uint64_t*)xcalloc(n_seqs, 8); isizes = (uint64_t*)calloc(n_seqs, 8);
for (i = 0, tot = 0; i != n_seqs; ++i) { for (i = 0, tot = 0; i != n_seqs; ++i) {
bwa_seq_t *p[2]; bwa_seq_t *p[2];
p[0] = seqs[0] + i; p[1] = seqs[1] + i; p[0] = seqs[0] + i; p[1] = seqs[1] + i;
@ -263,9 +267,9 @@ int bwa_cal_pac_pos_pe(const bntseq_t *bns, const char *prefix, bwt_t *const _bw
pe_data_t *d; pe_data_t *d;
aln_buf_t *buf[2]; aln_buf_t *buf[2];
d = (pe_data_t*)xcalloc(1, sizeof(pe_data_t)); d = (pe_data_t*)calloc(1, sizeof(pe_data_t));
buf[0] = (aln_buf_t*)xcalloc(n_seqs, sizeof(aln_buf_t)); buf[0] = (aln_buf_t*)calloc(n_seqs, sizeof(aln_buf_t));
buf[1] = (aln_buf_t*)xcalloc(n_seqs, sizeof(aln_buf_t)); buf[1] = (aln_buf_t*)calloc(n_seqs, sizeof(aln_buf_t));
if (_bwt == 0) { // load forward SA if (_bwt == 0) { // load forward SA
strcpy(str, prefix); strcat(str, ".bwt"); bwt = bwt_restore_bwt(str); strcpy(str, prefix); strcat(str, ".bwt"); bwt = bwt_restore_bwt(str);
@ -338,7 +342,7 @@ int bwa_cal_pac_pos_pe(const bntseq_t *bns, const char *prefix, bwt_t *const _bw
if (ret) { // not in the hash table; ret must equal 1 as we never remove elements if (ret) { // not in the hash table; ret must equal 1 as we never remove elements
poslist_t *z = &kh_val(g_hash, iter); poslist_t *z = &kh_val(g_hash, iter);
z->n = r->l - r->k + 1; z->n = r->l - r->k + 1;
z->a = (bwtint_t*)xmalloc(sizeof(bwtint_t) * z->n); z->a = (bwtint_t*)malloc(sizeof(bwtint_t) * z->n);
for (l = r->k; l <= r->l; ++l) { for (l = r->k; l <= r->l; ++l) {
int strand; int strand;
z->a[l - r->k] = bwa_sa2pos(bns, bwt, l, p[j]->len, &strand)<<1; z->a[l - r->k] = bwa_sa2pos(bns, bwt, l, p[j]->len, &strand)<<1;
@ -420,7 +424,7 @@ bwa_cigar_t *bwa_sw_core(bwtint_t l_pac, const ubyte_t *pacseq, int len, const u
if ((float)x/len >= 0.25 || len - x < SW_MIN_MATCH_LEN) return 0; if ((float)x/len >= 0.25 || len - x < SW_MIN_MATCH_LEN) return 0;
// get reference subsequence // get reference subsequence
ref_seq = (ubyte_t*)xcalloc(reglen, 1); ref_seq = (ubyte_t*)calloc(reglen, 1);
for (k = *beg, l = 0; l < reglen && k < l_pac; ++k) for (k = *beg, l = 0; l < reglen && k < l_pac; ++k)
ref_seq[l++] = pacseq[k>>2] >> ((~k&3)<<1) & 3; ref_seq[l++] = pacseq[k>>2] >> ((~k&3)<<1) & 3;
@ -453,7 +457,7 @@ bwa_cigar_t *bwa_sw_core(bwtint_t l_pac, const ubyte_t *pacseq, int len, const u
{ // update cigar and coordinate; { // update cigar and coordinate;
int start = r.qb, end = r.qe + 1; int start = r.qb, end = r.qe + 1;
*beg += r.tb; *beg += r.tb;
cigar = (bwa_cigar_t*)xrealloc(cigar, sizeof(bwa_cigar_t) * (*n_cigar + 2)); cigar = (bwa_cigar_t*)realloc(cigar, sizeof(bwa_cigar_t) * (*n_cigar + 2));
if (start) { if (start) {
memmove(cigar + 1, cigar, sizeof(bwa_cigar_t) * (*n_cigar)); memmove(cigar + 1, cigar, sizeof(bwa_cigar_t) * (*n_cigar));
cigar[0] = __cigar_create(3, start); cigar[0] = __cigar_create(3, start);
@ -497,7 +501,7 @@ ubyte_t *bwa_paired_sw(const bntseq_t *bns, const ubyte_t *_pacseq, int n_seqs,
// load reference sequence // load reference sequence
if (_pacseq == 0) { if (_pacseq == 0) {
pacseq = (ubyte_t*)xcalloc(bns->l_pac/4+1, 1); pacseq = (ubyte_t*)calloc(bns->l_pac/4+1, 1);
err_rewind(bns->fp_pac); err_rewind(bns->fp_pac);
err_fread_noeof(pacseq, 1, bns->l_pac/4+1, bns->fp_pac); err_fread_noeof(pacseq, 1, bns->l_pac/4+1, bns->fp_pac);
} else pacseq = (ubyte_t*)_pacseq; } else pacseq = (ubyte_t*)_pacseq;
@ -653,7 +657,7 @@ void bwa_sai2sam_pe_core(const char *prefix, char *const fn_sa[2], char *const f
if (popt->is_preload) { if (popt->is_preload) {
strcpy(str, prefix); strcat(str, ".bwt"); bwt = bwt_restore_bwt(str); strcpy(str, prefix); strcat(str, ".bwt"); bwt = bwt_restore_bwt(str);
strcpy(str, prefix); strcat(str, ".sa"); bwt_restore_sa(str, bwt); strcpy(str, prefix); strcat(str, ".sa"); bwt_restore_sa(str, bwt);
pac = (ubyte_t*)xcalloc(bns->l_pac/4+1, 1); pac = (ubyte_t*)calloc(bns->l_pac/4+1, 1);
err_rewind(bns->fp_pac); err_rewind(bns->fp_pac);
err_fread_noeof(pac, 1, bns->l_pac/4+1, bns->fp_pac); err_fread_noeof(pac, 1, bns->l_pac/4+1, bns->fp_pac);
} }

26
bwase.c
View File

@ -12,6 +12,10 @@
#include "bwa.h" #include "bwa.h"
#include "ksw.h" #include "ksw.h"
#ifdef USE_MALLOC_WRAPPERS
# include "malloc_wrap.h"
#endif
int g_log_n[256]; int g_log_n[256];
void bwa_print_sam_PG(); void bwa_print_sam_PG();
@ -59,7 +63,7 @@ void bwa_aln2seq_core(int n_aln, const bwt_aln1_t *aln, bwa_seq_t *s, int set_ma
* simply output all hits, but the following samples "rest" * simply output all hits, but the following samples "rest"
* number of random hits. */ * number of random hits. */
rest = n_occ > n_multi + 1? n_multi + 1 : n_occ; // find one additional for ->sa rest = n_occ > n_multi + 1? n_multi + 1 : n_occ; // find one additional for ->sa
s->multi = xcalloc(rest, sizeof(bwt_multi1_t)); s->multi = calloc(rest, sizeof(bwt_multi1_t));
for (k = 0; k < n_aln; ++k) { for (k = 0; k < n_aln; ++k) {
const bwt_aln1_t *q = aln + k; const bwt_aln1_t *q = aln + k;
if (q->l - q->k + 1 <= rest) { if (q->l - q->k + 1 <= rest) {
@ -184,7 +188,7 @@ bwa_cigar_t *bwa_refine_gapped_core(bwtint_t l_pac, const ubyte_t *pacseq, int l
if (rlen == 0) goto refine_gapped_err; if (rlen == 0) goto refine_gapped_err;
ksw_global(qle, &seq[len-qle], rlen, rseq, 5, mat, 5, 1, SW_BW, n_cigar, &cigar32); ksw_global(qle, &seq[len-qle], rlen, rseq, 5, mat, 5, 1, SW_BW, n_cigar, &cigar32);
if (qle < len) { // write soft clip if (qle < len) { // write soft clip
cigar = xrealloc(cigar, (*n_cigar + 1) * 4); cigar = realloc(cigar, (*n_cigar + 1) * 4);
memmove(cigar + 1, cigar, *n_cigar * 4); memmove(cigar + 1, cigar, *n_cigar * 4);
cigar[0] = (len - qle)<<4 | FROM_S; cigar[0] = (len - qle)<<4 | FROM_S;
++(*n_cigar); ++(*n_cigar);
@ -199,7 +203,7 @@ bwa_cigar_t *bwa_refine_gapped_core(bwtint_t l_pac, const ubyte_t *pacseq, int l
if (rlen == 0) goto refine_gapped_err; if (rlen == 0) goto refine_gapped_err;
ksw_global(qle, seq, rlen, rseq, 5, mat, 5, 1, SW_BW, n_cigar, &cigar32); // right extension ksw_global(qle, seq, rlen, rseq, 5, mat, 5, 1, SW_BW, n_cigar, &cigar32); // right extension
if (qle < len) { if (qle < len) {
cigar = xrealloc(cigar, (*n_cigar + 1) * 4); cigar = realloc(cigar, (*n_cigar + 1) * 4);
cigar[*n_cigar - 1] = (len - qle)<<4 | FROM_S; cigar[*n_cigar - 1] = (len - qle)<<4 | FROM_S;
++(*n_cigar); ++(*n_cigar);
} }
@ -265,7 +269,7 @@ char *bwa_cal_md1(int n_cigar, bwa_cigar_t *cigar, int len, bwtint_t pos, ubyte_
} }
ksprintf(str, "%d", u); ksprintf(str, "%d", u);
*_nm = nm; *_nm = nm;
return xstrdup(str->s); return strdup(str->s);
} }
void bwa_correct_trimmed(bwa_seq_t *s) void bwa_correct_trimmed(bwa_seq_t *s)
@ -277,11 +281,11 @@ void bwa_correct_trimmed(bwa_seq_t *s)
} else { } else {
if (s->cigar == 0) { if (s->cigar == 0) {
s->n_cigar = 2; s->n_cigar = 2;
s->cigar = xcalloc(s->n_cigar, sizeof(bwa_cigar_t)); s->cigar = calloc(s->n_cigar, sizeof(bwa_cigar_t));
s->cigar[0] = __cigar_create(0, s->len); s->cigar[0] = __cigar_create(0, s->len);
} else { } else {
++s->n_cigar; ++s->n_cigar;
s->cigar = xrealloc(s->cigar, s->n_cigar * sizeof(bwa_cigar_t)); s->cigar = realloc(s->cigar, s->n_cigar * sizeof(bwa_cigar_t));
} }
s->cigar[s->n_cigar-1] = __cigar_create(3, (s->full_len - s->len)); s->cigar[s->n_cigar-1] = __cigar_create(3, (s->full_len - s->len));
} }
@ -291,11 +295,11 @@ void bwa_correct_trimmed(bwa_seq_t *s)
} else { } else {
if (s->cigar == 0) { if (s->cigar == 0) {
s->n_cigar = 2; s->n_cigar = 2;
s->cigar = xcalloc(s->n_cigar, sizeof(bwa_cigar_t)); s->cigar = calloc(s->n_cigar, sizeof(bwa_cigar_t));
s->cigar[1] = __cigar_create(0, s->len); s->cigar[1] = __cigar_create(0, s->len);
} else { } else {
++s->n_cigar; ++s->n_cigar;
s->cigar = xrealloc(s->cigar, s->n_cigar * sizeof(bwa_cigar_t)); s->cigar = realloc(s->cigar, s->n_cigar * sizeof(bwa_cigar_t));
memmove(s->cigar + 1, s->cigar, (s->n_cigar-1) * sizeof(bwa_cigar_t)); memmove(s->cigar + 1, s->cigar, (s->n_cigar-1) * sizeof(bwa_cigar_t));
} }
s->cigar[0] = __cigar_create(3, (s->full_len - s->len)); s->cigar[0] = __cigar_create(3, (s->full_len - s->len));
@ -311,7 +315,7 @@ void bwa_refine_gapped(const bntseq_t *bns, int n_seqs, bwa_seq_t *seqs, ubyte_t
kstring_t *str; kstring_t *str;
if (!_pacseq) { if (!_pacseq) {
pacseq = (ubyte_t*)xcalloc(bns->l_pac/4+1, 1); pacseq = (ubyte_t*)calloc(bns->l_pac/4+1, 1);
err_rewind(bns->fp_pac); err_rewind(bns->fp_pac);
err_fread_noeof(pacseq, 1, bns->l_pac/4+1, bns->fp_pac); err_fread_noeof(pacseq, 1, bns->l_pac/4+1, bns->fp_pac);
} else pacseq = _pacseq; } else pacseq = _pacseq;
@ -330,7 +334,7 @@ void bwa_refine_gapped(const bntseq_t *bns, int n_seqs, bwa_seq_t *seqs, ubyte_t
if (s->cigar == 0) s->type = BWA_TYPE_NO_MATCH; if (s->cigar == 0) s->type = BWA_TYPE_NO_MATCH;
} }
// generate MD tag // generate MD tag
str = (kstring_t*)xcalloc(1, sizeof(kstring_t)); str = (kstring_t*)calloc(1, sizeof(kstring_t));
for (i = 0; i != n_seqs; ++i) { for (i = 0; i != n_seqs; ++i) {
bwa_seq_t *s = seqs + i; bwa_seq_t *s = seqs + i;
if (s->type != BWA_TYPE_NO_MATCH) { if (s->type != BWA_TYPE_NO_MATCH) {
@ -559,7 +563,7 @@ void bwa_sai2sam_se_core(const char *prefix, const char *fn_sa, const char *fn_f
err_fread_noeof(&n_aln, 4, 1, fp_sa); err_fread_noeof(&n_aln, 4, 1, fp_sa);
if (n_aln > m_aln) { if (n_aln > m_aln) {
m_aln = n_aln; m_aln = n_aln;
aln = (bwt_aln1_t*)xrealloc(aln, sizeof(bwt_aln1_t) * m_aln); aln = (bwt_aln1_t*)realloc(aln, sizeof(bwt_aln1_t) * m_aln);
} }
err_fread_noeof(aln, sizeof(bwt_aln1_t), n_aln, fp_sa); err_fread_noeof(aln, sizeof(bwt_aln1_t), n_aln, fp_sa);
bwa_aln2seq_core(n_aln, aln, p, 1, n_occ); bwa_aln2seq_core(n_aln, aln, p, 1, n_occ);

28
bwaseqio.c
View File

@ -7,6 +7,10 @@
#include "kseq.h" #include "kseq.h"
KSEQ_DECLARE(gzFile) KSEQ_DECLARE(gzFile)
#ifdef USE_MALLOC_WRAPPERS
# include "malloc_wrap.h"
#endif
extern unsigned char nst_nt4_table[256]; extern unsigned char nst_nt4_table[256];
static char bam_nt16_nt4_table[] = { 4, 0, 1, 4, 2, 4, 4, 4, 3, 4, 4, 4, 4, 4, 4, 4 }; static char bam_nt16_nt4_table[] = { 4, 0, 1, 4, 2, 4, 4, 4, 3, 4, 4, 4, 4, 4, 4, 4 };
@ -22,7 +26,7 @@ bwa_seqio_t *bwa_bam_open(const char *fn, int which)
{ {
bwa_seqio_t *bs; bwa_seqio_t *bs;
bam_header_t *h; bam_header_t *h;
bs = (bwa_seqio_t*)xcalloc(1, sizeof(bwa_seqio_t)); bs = (bwa_seqio_t*)calloc(1, sizeof(bwa_seqio_t));
bs->is_bam = 1; bs->is_bam = 1;
bs->which = which; bs->which = which;
bs->fp = bam_open(fn, "r"); bs->fp = bam_open(fn, "r");
@ -35,7 +39,7 @@ bwa_seqio_t *bwa_seq_open(const char *fn)
{ {
gzFile fp; gzFile fp;
bwa_seqio_t *bs; bwa_seqio_t *bs;
bs = (bwa_seqio_t*)xcalloc(1, sizeof(bwa_seqio_t)); bs = (bwa_seqio_t*)calloc(1, sizeof(bwa_seqio_t));
fp = xzopen(fn, "r"); fp = xzopen(fn, "r");
bs->ks = kseq_init(fp); bs->ks = kseq_init(fp);
return bs; return bs;
@ -93,7 +97,7 @@ static bwa_seq_t *bwa_read_bam(bwa_seqio_t *bs, int n_needed, int *n, int is_com
b = bam_init1(); b = bam_init1();
n_seqs = 0; n_seqs = 0;
seqs = (bwa_seq_t*)xcalloc(n_needed, sizeof(bwa_seq_t)); seqs = (bwa_seq_t*)calloc(n_needed, sizeof(bwa_seq_t));
while (bam_read1(bs->fp, b) >= 0) { while (bam_read1(bs->fp, b) >= 0) {
uint8_t *s, *q; uint8_t *s, *q;
int go = 0; int go = 0;
@ -108,8 +112,8 @@ static bwa_seq_t *bwa_read_bam(bwa_seqio_t *bs, int n_needed, int *n, int is_com
p->full_len = p->clip_len = p->len = l; p->full_len = p->clip_len = p->len = l;
n_tot += p->full_len; n_tot += p->full_len;
s = bam1_seq(b); q = bam1_qual(b); s = bam1_seq(b); q = bam1_qual(b);
p->seq = (ubyte_t*)xcalloc(p->len + 1, 1); p->seq = (ubyte_t*)calloc(p->len + 1, 1);
p->qual = (ubyte_t*)xcalloc(p->len + 1, 1); p->qual = (ubyte_t*)calloc(p->len + 1, 1);
for (i = 0; i != p->full_len; ++i) { for (i = 0; i != p->full_len; ++i) {
p->seq[i] = bam_nt16_nt4_table[(int)bam1_seqi(s, i)]; p->seq[i] = bam_nt16_nt4_table[(int)bam1_seqi(s, i)];
p->qual[i] = q[i] + 33 < 126? q[i] + 33 : 126; p->qual[i] = q[i] + 33 < 126? q[i] + 33 : 126;
@ -119,11 +123,11 @@ static bwa_seq_t *bwa_read_bam(bwa_seqio_t *bs, int n_needed, int *n, int is_com
seq_reverse(p->len, p->qual, 0); seq_reverse(p->len, p->qual, 0);
} }
if (trim_qual >= 1) n_trimmed += bwa_trim_read(trim_qual, p); if (trim_qual >= 1) n_trimmed += bwa_trim_read(trim_qual, p);
p->rseq = (ubyte_t*)xcalloc(p->full_len, 1); p->rseq = (ubyte_t*)calloc(p->full_len, 1);
memcpy(p->rseq, p->seq, p->len); memcpy(p->rseq, p->seq, p->len);
seq_reverse(p->len, p->seq, 0); // *IMPORTANT*: will be reversed back in bwa_refine_gapped() seq_reverse(p->len, p->seq, 0); // *IMPORTANT*: will be reversed back in bwa_refine_gapped()
seq_reverse(p->len, p->rseq, is_comp); seq_reverse(p->len, p->rseq, is_comp);
p->name = xstrdup((const char*)bam1_qname(b)); p->name = strdup((const char*)bam1_qname(b));
if (n_seqs == n_needed) break; if (n_seqs == n_needed) break;
} }
*n = n_seqs; *n = n_seqs;
@ -153,7 +157,7 @@ bwa_seq_t *bwa_read_seq(bwa_seqio_t *bs, int n_needed, int *n, int mode, int tri
} }
if (bs->is_bam) return bwa_read_bam(bs, n_needed, n, is_comp, trim_qual); // l_bc has no effect for BAM input if (bs->is_bam) return bwa_read_bam(bs, n_needed, n, is_comp, trim_qual); // l_bc has no effect for BAM input
n_seqs = 0; n_seqs = 0;
seqs = (bwa_seq_t*)xcalloc(n_needed, sizeof(bwa_seq_t)); seqs = (bwa_seq_t*)calloc(n_needed, sizeof(bwa_seq_t));
while ((l = kseq_read(seq)) >= 0) { while ((l = kseq_read(seq)) >= 0) {
if ((mode & BWA_MODE_CFY) && (seq->comment.l != 0)) { if ((mode & BWA_MODE_CFY) && (seq->comment.l != 0)) {
// skip reads that are marked to be filtered by Casava // skip reads that are marked to be filtered by Casava
@ -184,18 +188,18 @@ bwa_seq_t *bwa_read_seq(bwa_seqio_t *bs, int n_needed, int *n, int mode, int tri
p->qual = 0; p->qual = 0;
p->full_len = p->clip_len = p->len = l; p->full_len = p->clip_len = p->len = l;
n_tot += p->full_len; n_tot += p->full_len;
p->seq = (ubyte_t*)xcalloc(p->full_len, 1); p->seq = (ubyte_t*)calloc(p->full_len, 1);
for (i = 0; i != p->full_len; ++i) for (i = 0; i != p->full_len; ++i)
p->seq[i] = nst_nt4_table[(int)seq->seq.s[i]]; p->seq[i] = nst_nt4_table[(int)seq->seq.s[i]];
if (seq->qual.l) { // copy quality if (seq->qual.l) { // copy quality
p->qual = (ubyte_t*)xstrdup((char*)seq->qual.s); p->qual = (ubyte_t*)strdup((char*)seq->qual.s);
if (trim_qual >= 1) n_trimmed += bwa_trim_read(trim_qual, p); if (trim_qual >= 1) n_trimmed += bwa_trim_read(trim_qual, p);
} }
p->rseq = (ubyte_t*)xcalloc(p->full_len, 1); p->rseq = (ubyte_t*)calloc(p->full_len, 1);
memcpy(p->rseq, p->seq, p->len); memcpy(p->rseq, p->seq, p->len);
seq_reverse(p->len, p->seq, 0); // *IMPORTANT*: will be reversed back in bwa_refine_gapped() seq_reverse(p->len, p->seq, 0); // *IMPORTANT*: will be reversed back in bwa_refine_gapped()
seq_reverse(p->len, p->rseq, is_comp); seq_reverse(p->len, p->rseq, is_comp);
p->name = xstrdup((const char*)seq->name.s); p->name = strdup((const char*)seq->name.s);
{ // trim /[12]$ { // trim /[12]$
int t = strlen(p->name); int t = strlen(p->name);
if (t > 2 && p->name[t-2] == '/' && (p->name[t-1] == '1' || p->name[t-1] == '2')) p->name[t-2] = '\0'; if (t > 2 && p->name[t-2] == '/' && (p->name[t-1] == '1' || p->name[t-1] == '2')) p->name[t-2] = '\0';

12
bwt.c
View File

@ -34,6 +34,10 @@
#include "bwt.h" #include "bwt.h"
#include "kvec.h" #include "kvec.h"
#ifdef USE_MALLOC_WRAPPERS
# include "malloc_wrap.h"
#endif
void bwt_gen_cnt_table(bwt_t *bwt) void bwt_gen_cnt_table(bwt_t *bwt)
{ {
int i, j; int i, j;
@ -66,7 +70,7 @@ void bwt_cal_sa(bwt_t *bwt, int intv)
if (bwt->sa) free(bwt->sa); if (bwt->sa) free(bwt->sa);
bwt->sa_intv = intv; bwt->sa_intv = intv;
bwt->n_sa = (bwt->seq_len + intv) / intv; bwt->n_sa = (bwt->seq_len + intv) / intv;
bwt->sa = (bwtint_t*)xcalloc(bwt->n_sa, sizeof(bwtint_t)); bwt->sa = (bwtint_t*)calloc(bwt->n_sa, sizeof(bwtint_t));
// calculate SA value // calculate SA value
isa = 0; sa = bwt->seq_len; isa = 0; sa = bwt->seq_len;
for (i = 0; i < bwt->seq_len; ++i) { for (i = 0; i < bwt->seq_len; ++i) {
@ -397,7 +401,7 @@ void bwt_restore_sa(const char *fn, bwt_t *bwt)
xassert(primary == bwt->seq_len, "SA-BWT inconsistency: seq_len is not the same."); xassert(primary == bwt->seq_len, "SA-BWT inconsistency: seq_len is not the same.");
bwt->n_sa = (bwt->seq_len + bwt->sa_intv) / bwt->sa_intv; bwt->n_sa = (bwt->seq_len + bwt->sa_intv) / bwt->sa_intv;
bwt->sa = (bwtint_t*)xcalloc(bwt->n_sa, sizeof(bwtint_t)); bwt->sa = (bwtint_t*)calloc(bwt->n_sa, sizeof(bwtint_t));
bwt->sa[0] = -1; bwt->sa[0] = -1;
fread_fix(fp, sizeof(bwtint_t) * (bwt->n_sa - 1), bwt->sa + 1); fread_fix(fp, sizeof(bwtint_t) * (bwt->n_sa - 1), bwt->sa + 1);
@ -409,11 +413,11 @@ bwt_t *bwt_restore_bwt(const char *fn)
bwt_t *bwt; bwt_t *bwt;
FILE *fp; FILE *fp;
bwt = (bwt_t*)xcalloc(1, sizeof(bwt_t)); bwt = (bwt_t*)calloc(1, sizeof(bwt_t));
fp = xopen(fn, "rb"); fp = xopen(fn, "rb");
err_fseek(fp, 0, SEEK_END); err_fseek(fp, 0, SEEK_END);
bwt->bwt_size = (err_ftell(fp) - sizeof(bwtint_t) * 5) >> 2; bwt->bwt_size = (err_ftell(fp) - sizeof(bwtint_t) * 5) >> 2;
bwt->bwt = (uint32_t*)xcalloc(bwt->bwt_size, 4); bwt->bwt = (uint32_t*)calloc(bwt->bwt_size, 4);
err_fseek(fp, 0, SEEK_SET); err_fseek(fp, 0, SEEK_SET);
err_fread_noeof(&bwt->primary, sizeof(bwtint_t), 1, fp); err_fread_noeof(&bwt->primary, sizeof(bwtint_t), 1, fp);
err_fread_noeof(bwt->L2+1, sizeof(bwtint_t), 4, fp); err_fread_noeof(bwt->L2+1, sizeof(bwtint_t), 4, fp);

20
bwt_gen.c
View File

@ -30,6 +30,10 @@
#include "QSufSort.h" #include "QSufSort.h"
#include "utils.h" #include "utils.h"
#ifdef USE_MALLOC_WRAPPERS
# include "malloc_wrap.h"
#endif
typedef uint64_t bgint_t; typedef uint64_t bgint_t;
typedef int64_t sbgint_t; typedef int64_t sbgint_t;
@ -320,25 +324,25 @@ BWT *BWTCreate(const bgint_t textLength, unsigned int *decodeTable)
{ {
BWT *bwt; BWT *bwt;
bwt = (BWT*)xcalloc(1, sizeof(BWT)); bwt = (BWT*)calloc(1, sizeof(BWT));
bwt->textLength = 0; bwt->textLength = 0;
bwt->cumulativeFreq = (bgint_t*)xcalloc((ALPHABET_SIZE + 1), sizeof(bgint_t)); bwt->cumulativeFreq = (bgint_t*)calloc((ALPHABET_SIZE + 1), sizeof(bgint_t));
initializeVAL_bg(bwt->cumulativeFreq, ALPHABET_SIZE + 1, 0); initializeVAL_bg(bwt->cumulativeFreq, ALPHABET_SIZE + 1, 0);
bwt->bwtSizeInWord = 0; bwt->bwtSizeInWord = 0;
// Generate decode tables // Generate decode tables
if (decodeTable == NULL) { if (decodeTable == NULL) {
bwt->decodeTable = (unsigned*)xcalloc(DNA_OCC_CNT_TABLE_SIZE_IN_WORD, sizeof(unsigned int)); bwt->decodeTable = (unsigned*)calloc(DNA_OCC_CNT_TABLE_SIZE_IN_WORD, sizeof(unsigned int));
GenerateDNAOccCountTable(bwt->decodeTable); GenerateDNAOccCountTable(bwt->decodeTable);
} else { } else {
bwt->decodeTable = decodeTable; bwt->decodeTable = decodeTable;
} }
bwt->occMajorSizeInWord = BWTOccValueMajorSizeInWord(textLength); bwt->occMajorSizeInWord = BWTOccValueMajorSizeInWord(textLength);
bwt->occValueMajor = (bgint_t*)xcalloc(bwt->occMajorSizeInWord, sizeof(bgint_t)); bwt->occValueMajor = (bgint_t*)calloc(bwt->occMajorSizeInWord, sizeof(bgint_t));
bwt->occSizeInWord = 0; bwt->occSizeInWord = 0;
bwt->occValue = NULL; bwt->occValue = NULL;
@ -354,16 +358,16 @@ BWTInc *BWTIncCreate(const bgint_t textLength, unsigned int initialMaxBuildSize,
if (textLength < incMaxBuildSize) incMaxBuildSize = textLength; if (textLength < incMaxBuildSize) incMaxBuildSize = textLength;
if (textLength < initialMaxBuildSize) initialMaxBuildSize = textLength; if (textLength < initialMaxBuildSize) initialMaxBuildSize = textLength;
bwtInc = (BWTInc*)xcalloc(1, sizeof(BWTInc)); bwtInc = (BWTInc*)calloc(1, sizeof(BWTInc));
bwtInc->numberOfIterationDone = 0; bwtInc->numberOfIterationDone = 0;
bwtInc->bwt = BWTCreate(textLength, NULL); bwtInc->bwt = BWTCreate(textLength, NULL);
bwtInc->initialMaxBuildSize = initialMaxBuildSize; bwtInc->initialMaxBuildSize = initialMaxBuildSize;
bwtInc->incMaxBuildSize = incMaxBuildSize; bwtInc->incMaxBuildSize = incMaxBuildSize;
bwtInc->cumulativeCountInCurrentBuild = (bgint_t*)xcalloc((ALPHABET_SIZE + 1), sizeof(bgint_t)); bwtInc->cumulativeCountInCurrentBuild = (bgint_t*)calloc((ALPHABET_SIZE + 1), sizeof(bgint_t));
initializeVAL_bg(bwtInc->cumulativeCountInCurrentBuild, ALPHABET_SIZE + 1, 0); initializeVAL_bg(bwtInc->cumulativeCountInCurrentBuild, ALPHABET_SIZE + 1, 0);
// Build frequently accessed data // Build frequently accessed data
bwtInc->packedShift = (unsigned*)xcalloc(CHAR_PER_WORD, sizeof(unsigned int)); bwtInc->packedShift = (unsigned*)calloc(CHAR_PER_WORD, sizeof(unsigned int));
for (i=0; i<CHAR_PER_WORD; i++) for (i=0; i<CHAR_PER_WORD; i++)
bwtInc->packedShift[i] = BITS_IN_WORD - (i+1) * BIT_PER_CHAR; bwtInc->packedShift[i] = BITS_IN_WORD - (i+1) * BIT_PER_CHAR;
@ -373,7 +377,7 @@ BWTInc *BWTIncCreate(const bgint_t textLength, unsigned int initialMaxBuildSize,
+ incMaxBuildSize/5 * 3 * (sizeof(bgint_t) / 4); // space for the 3 temporary arrays in each iteration + incMaxBuildSize/5 * 3 * (sizeof(bgint_t) / 4); // space for the 3 temporary arrays in each iteration
if (bwtInc->availableWord < MIN_AVAILABLE_WORD) bwtInc->availableWord = MIN_AVAILABLE_WORD; // lh3: otherwise segfaul when availableWord is too small if (bwtInc->availableWord < MIN_AVAILABLE_WORD) bwtInc->availableWord = MIN_AVAILABLE_WORD; // lh3: otherwise segfaul when availableWord is too small
fprintf(stderr, "[%s] textLength=%ld, availableWord=%ld\n", __func__, (long)textLength, (long)bwtInc->availableWord); fprintf(stderr, "[%s] textLength=%ld, availableWord=%ld\n", __func__, (long)textLength, (long)bwtInc->availableWord);
bwtInc->workingMemory = (unsigned*)xcalloc(bwtInc->availableWord, BYTES_IN_WORD); bwtInc->workingMemory = (unsigned*)calloc(bwtInc->availableWord, BYTES_IN_WORD);
return bwtInc; return bwtInc;
} }

15
bwt_lite.c
View File

@ -2,7 +2,10 @@
#include <string.h> #include <string.h>
#include <stdio.h> #include <stdio.h>
#include "bwt_lite.h" #include "bwt_lite.h"
#include "utils.h"
#ifdef USE_MALLOC_WRAPPERS
# include "malloc_wrap.h"
#endif
int is_sa(const uint8_t *T, uint32_t *SA, int n); int is_sa(const uint8_t *T, uint32_t *SA, int n);
int is_bwt(uint8_t *T, int n); int is_bwt(uint8_t *T, int n);
@ -11,21 +14,21 @@ bwtl_t *bwtl_seq2bwtl(int len, const uint8_t *seq)
{ {
bwtl_t *b; bwtl_t *b;
int i; int i;
b = (bwtl_t*)xcalloc(1, sizeof(bwtl_t)); b = (bwtl_t*)calloc(1, sizeof(bwtl_t));
b->seq_len = len; b->seq_len = len;
{ // calculate b->bwt { // calculate b->bwt
uint8_t *s; uint8_t *s;
b->sa = (uint32_t*)xcalloc(len + 1, 4); b->sa = (uint32_t*)calloc(len + 1, 4);
is_sa(seq, b->sa, len); is_sa(seq, b->sa, len);
s = (uint8_t*)xcalloc(len + 1, 1); s = (uint8_t*)calloc(len + 1, 1);
for (i = 0; i <= len; ++i) { for (i = 0; i <= len; ++i) {
if (b->sa[i] == 0) b->primary = i; if (b->sa[i] == 0) b->primary = i;
else s[i] = seq[b->sa[i] - 1]; else s[i] = seq[b->sa[i] - 1];
} }
for (i = b->primary; i < len; ++i) s[i] = s[i + 1]; for (i = b->primary; i < len; ++i) s[i] = s[i + 1];
b->bwt_size = (len + 15) / 16; b->bwt_size = (len + 15) / 16;
b->bwt = (uint32_t*)xcalloc(b->bwt_size, 4); b->bwt = (uint32_t*)calloc(b->bwt_size, 4);
for (i = 0; i < len; ++i) for (i = 0; i < len; ++i)
b->bwt[i>>4] |= s[i] << ((15 - (i&15)) << 1); b->bwt[i>>4] |= s[i] << ((15 - (i&15)) << 1);
free(s); free(s);
@ -33,7 +36,7 @@ bwtl_t *bwtl_seq2bwtl(int len, const uint8_t *seq)
{ // calculate b->occ { // calculate b->occ
uint32_t c[4]; uint32_t c[4];
b->n_occ = (len + 15) / 16 * 4; b->n_occ = (len + 15) / 16 * 4;
b->occ = (uint32_t*)xcalloc(b->n_occ, 4); b->occ = (uint32_t*)calloc(b->n_occ, 4);
memset(c, 0, 16); memset(c, 0, 16);
for (i = 0; i < len; ++i) { for (i = 0; i < len; ++i) {
if (i % 16 == 0) if (i % 16 == 0)

16
bwtaln.c
View File

@ -17,10 +17,14 @@
#include <pthread.h> #include <pthread.h>
#endif #endif
#ifdef USE_MALLOC_WRAPPERS
# include "malloc_wrap.h"
#endif
gap_opt_t *gap_init_opt() gap_opt_t *gap_init_opt()
{ {
gap_opt_t *o; gap_opt_t *o;
o = (gap_opt_t*)xcalloc(1, sizeof(gap_opt_t)); o = (gap_opt_t*)calloc(1, sizeof(gap_opt_t));
/* IMPORTANT: s_mm*10 should be about the average base error /* IMPORTANT: s_mm*10 should be about the average base error
rate. Voilating this requirement will break pairing! */ rate. Voilating this requirement will break pairing! */
o->s_mm = 3; o->s_gapo = 11; o->s_gape = 4; o->s_mm = 3; o->s_gapo = 11; o->s_gape = 4;
@ -90,7 +94,7 @@ void bwa_cal_sa_reg_gap(int tid, bwt_t *const bwt, int n_seqs, bwa_seq_t *seqs,
if (local_opt.max_diff < local_opt.max_gapo) local_opt.max_gapo = local_opt.max_diff; if (local_opt.max_diff < local_opt.max_gapo) local_opt.max_gapo = local_opt.max_diff;
stack = gap_init_stack(local_opt.max_diff, local_opt.max_gapo, local_opt.max_gape, &local_opt); stack = gap_init_stack(local_opt.max_diff, local_opt.max_gapo, local_opt.max_gape, &local_opt);
seed_w = (bwt_width_t*)xcalloc(opt->seed_len+1, sizeof(bwt_width_t)); seed_w = (bwt_width_t*)calloc(opt->seed_len+1, sizeof(bwt_width_t));
w = 0; w = 0;
for (i = 0; i != n_seqs; ++i) { for (i = 0; i != n_seqs; ++i) {
bwa_seq_t *p = seqs + i; bwa_seq_t *p = seqs + i;
@ -100,7 +104,7 @@ void bwa_cal_sa_reg_gap(int tid, bwt_t *const bwt, int n_seqs, bwa_seq_t *seqs,
p->sa = 0; p->type = BWA_TYPE_NO_MATCH; p->c1 = p->c2 = 0; p->n_aln = 0; p->aln = 0; p->sa = 0; p->type = BWA_TYPE_NO_MATCH; p->c1 = p->c2 = 0; p->n_aln = 0; p->aln = 0;
if (max_l < p->len) { if (max_l < p->len) {
max_l = p->len; max_l = p->len;
w = (bwt_width_t*)xrealloc(w, (max_l + 1) * sizeof(bwt_width_t)); w = (bwt_width_t*)realloc(w, (max_l + 1) * sizeof(bwt_width_t));
memset(w, 0, (max_l + 1) * sizeof(bwt_width_t)); memset(w, 0, (max_l + 1) * sizeof(bwt_width_t));
} }
bwt_cal_width(bwt, p->len, p->seq, w); bwt_cal_width(bwt, p->len, p->seq, w);
@ -163,7 +167,7 @@ void bwa_aln_core(const char *prefix, const char *fn_fa, const gap_opt_t *opt)
ks = bwa_open_reads(opt->mode, fn_fa); ks = bwa_open_reads(opt->mode, fn_fa);
{ // load BWT { // load BWT
char *str = (char*)xcalloc(strlen(prefix) + 10, 1); char *str = (char*)calloc(strlen(prefix) + 10, 1);
strcpy(str, prefix); strcat(str, ".bwt"); bwt = bwt_restore_bwt(str); strcpy(str, prefix); strcat(str, ".bwt"); bwt = bwt_restore_bwt(str);
free(str); free(str);
} }
@ -186,8 +190,8 @@ void bwa_aln_core(const char *prefix, const char *fn_fa, const gap_opt_t *opt)
int j; int j;
pthread_attr_init(&attr); pthread_attr_init(&attr);
pthread_attr_setdetachstate(&attr, PTHREAD_CREATE_JOINABLE); pthread_attr_setdetachstate(&attr, PTHREAD_CREATE_JOINABLE);
data = (thread_aux_t*)xcalloc(opt->n_threads, sizeof(thread_aux_t)); data = (thread_aux_t*)calloc(opt->n_threads, sizeof(thread_aux_t));
tid = (pthread_t*)xcalloc(opt->n_threads, sizeof(pthread_t)); tid = (pthread_t*)calloc(opt->n_threads, sizeof(pthread_t));
for (j = 0; j < opt->n_threads; ++j) { for (j = 0; j < opt->n_threads; ++j) {
data[j].tid = j; data[j].bwt = bwt; data[j].tid = j; data[j].bwt = bwt;
data[j].n_seqs = n_seqs; data[j].seqs = seqs; data[j].opt = opt; data[j].n_seqs = n_seqs; data[j].seqs = seqs; data[j].opt = opt;

15
bwtgap.c
View File

@ -3,7 +3,10 @@
#include <string.h> #include <string.h>
#include "bwtgap.h" #include "bwtgap.h"
#include "bwtaln.h" #include "bwtaln.h"
#include "utils.h"
#ifdef USE_MALLOC_WRAPPERS
# include "malloc_wrap.h"
#endif
#define STATE_M 0 #define STATE_M 0
#define STATE_I 1 #define STATE_I 1
@ -14,9 +17,9 @@
gap_stack_t *gap_init_stack2(int max_score) gap_stack_t *gap_init_stack2(int max_score)
{ {
gap_stack_t *stack; gap_stack_t *stack;
stack = (gap_stack_t*)xcalloc(1, sizeof(gap_stack_t)); stack = (gap_stack_t*)calloc(1, sizeof(gap_stack_t));
stack->n_stacks = max_score; stack->n_stacks = max_score;
stack->stacks = (gap_stack1_t*)xcalloc(stack->n_stacks, sizeof(gap_stack1_t)); stack->stacks = (gap_stack1_t*)calloc(stack->n_stacks, sizeof(gap_stack1_t));
return stack; return stack;
} }
@ -52,7 +55,7 @@ static inline void gap_push(gap_stack_t *stack, int i, bwtint_t k, bwtint_t l, i
q = stack->stacks + score; q = stack->stacks + score;
if (q->n_entries == q->m_entries) { if (q->n_entries == q->m_entries) {
q->m_entries = q->m_entries? q->m_entries<<1 : 4; q->m_entries = q->m_entries? q->m_entries<<1 : 4;
q->stack = (gap_entry_t*)xrealloc(q->stack, sizeof(gap_entry_t) * q->m_entries); q->stack = (gap_entry_t*)realloc(q->stack, sizeof(gap_entry_t) * q->m_entries);
} }
p = q->stack + q->n_entries; p = q->stack + q->n_entries;
p->info = (u_int32_t)score<<21 | i; p->k = k; p->l = l; p->info = (u_int32_t)score<<21 | i; p->k = k; p->l = l;
@ -111,7 +114,7 @@ bwt_aln1_t *bwt_match_gap(bwt_t *const bwt, int len, const ubyte_t *seq, bwt_wid
bwt_aln1_t *aln; bwt_aln1_t *aln;
m_aln = 4; n_aln = 0; m_aln = 4; n_aln = 0;
aln = (bwt_aln1_t*)xcalloc(m_aln, sizeof(bwt_aln1_t)); aln = (bwt_aln1_t*)calloc(m_aln, sizeof(bwt_aln1_t));
// check whether there are too many N // check whether there are too many N
for (j = _j = 0; j < len; ++j) for (j = _j = 0; j < len; ++j)
@ -178,7 +181,7 @@ bwt_aln1_t *bwt_match_gap(bwt_t *const bwt, int len, const ubyte_t *seq, bwt_wid
gap_shadow(l - k + 1, len, bwt->seq_len, e.last_diff_pos, width); gap_shadow(l - k + 1, len, bwt->seq_len, e.last_diff_pos, width);
if (n_aln == m_aln) { if (n_aln == m_aln) {
m_aln <<= 1; m_aln <<= 1;
aln = (bwt_aln1_t*)xrealloc(aln, m_aln * sizeof(bwt_aln1_t)); aln = (bwt_aln1_t*)realloc(aln, m_aln * sizeof(bwt_aln1_t));
memset(aln + m_aln/2, 0, m_aln/2*sizeof(bwt_aln1_t)); memset(aln + m_aln/2, 0, m_aln/2*sizeof(bwt_aln1_t));
} }
p = aln + n_aln; p = aln + n_aln;

25
bwtindex.c
View File

@ -39,6 +39,11 @@
#include "divsufsort.h" #include "divsufsort.h"
#endif #endif
#ifdef USE_MALLOC_WRAPPERS
# include "malloc_wrap.h"
#endif
int is_bwt(ubyte_t *T, int n); int is_bwt(ubyte_t *T, int n);
int64_t bwa_seq_len(const char *fn_pac) int64_t bwa_seq_len(const char *fn_pac)
@ -62,18 +67,18 @@ bwt_t *bwt_pac2bwt(const char *fn_pac, int use_is)
FILE *fp; FILE *fp;
// initialization // initialization
bwt = (bwt_t*)xcalloc(1, sizeof(bwt_t)); bwt = (bwt_t*)calloc(1, sizeof(bwt_t));
bwt->seq_len = bwa_seq_len(fn_pac); bwt->seq_len = bwa_seq_len(fn_pac);
bwt->bwt_size = (bwt->seq_len + 15) >> 4; bwt->bwt_size = (bwt->seq_len + 15) >> 4;
fp = xopen(fn_pac, "rb"); fp = xopen(fn_pac, "rb");
// prepare sequence // prepare sequence
pac_size = (bwt->seq_len>>2) + ((bwt->seq_len&3) == 0? 0 : 1); pac_size = (bwt->seq_len>>2) + ((bwt->seq_len&3) == 0? 0 : 1);
buf2 = (ubyte_t*)xcalloc(pac_size, 1); buf2 = (ubyte_t*)calloc(pac_size, 1);
err_fread_noeof(buf2, 1, pac_size, fp); err_fread_noeof(buf2, 1, pac_size, fp);
err_fclose(fp); err_fclose(fp);
memset(bwt->L2, 0, 5 * 4); memset(bwt->L2, 0, 5 * 4);
buf = (ubyte_t*)xcalloc(bwt->seq_len + 1, 1); buf = (ubyte_t*)calloc(bwt->seq_len + 1, 1);
for (i = 0; i < bwt->seq_len; ++i) { for (i = 0; i < bwt->seq_len; ++i) {
buf[i] = buf2[i>>2] >> ((3 - (i&3)) << 1) & 3; buf[i] = buf2[i>>2] >> ((3 - (i&3)) << 1) & 3;
++bwt->L2[1+buf[i]]; ++bwt->L2[1+buf[i]];
@ -91,7 +96,7 @@ bwt_t *bwt_pac2bwt(const char *fn_pac, int use_is)
err_fatal_simple("libdivsufsort is not compiled in."); err_fatal_simple("libdivsufsort is not compiled in.");
#endif #endif
} }
bwt->bwt = (u_int32_t*)xcalloc(bwt->bwt_size, 4); bwt->bwt = (u_int32_t*)calloc(bwt->bwt_size, 4);
for (i = 0; i < bwt->seq_len; ++i) for (i = 0; i < bwt->seq_len; ++i)
bwt->bwt[i>>4] |= buf[i] << ((15 - (i&15)) << 1); bwt->bwt[i>>4] |= buf[i] << ((15 - (i&15)) << 1);
free(buf); free(buf);
@ -127,7 +132,7 @@ void bwt_bwtupdate_core(bwt_t *bwt)
n_occ = (bwt->seq_len + OCC_INTERVAL - 1) / OCC_INTERVAL + 1; n_occ = (bwt->seq_len + OCC_INTERVAL - 1) / OCC_INTERVAL + 1;
bwt->bwt_size += n_occ * sizeof(bwtint_t); // the new size bwt->bwt_size += n_occ * sizeof(bwtint_t); // the new size
buf = (uint32_t*)xcalloc(bwt->bwt_size, 4); // will be the new bwt buf = (uint32_t*)calloc(bwt->bwt_size, 4); // will be the new bwt
c[0] = c[1] = c[2] = c[3] = 0; c[0] = c[1] = c[2] = c[3] = 0;
for (i = k = 0; i < bwt->seq_len; ++i) { for (i = k = 0; i < bwt->seq_len; ++i) {
if (i % OCC_INTERVAL == 0) { if (i % OCC_INTERVAL == 0) {
@ -196,7 +201,7 @@ int bwa_index(int argc, char *argv[]) // the "index" command
else if (strcmp(optarg, "is") == 0) algo_type = 3; else if (strcmp(optarg, "is") == 0) algo_type = 3;
else err_fatal(__func__, "unknown algorithm: '%s'.", optarg); else err_fatal(__func__, "unknown algorithm: '%s'.", optarg);
break; break;
case 'p': prefix = xstrdup(optarg); break; case 'p': prefix = strdup(optarg); break;
case '6': is_64 = 1; break; case '6': is_64 = 1; break;
default: return 1; default: return 1;
} }
@ -215,13 +220,13 @@ int bwa_index(int argc, char *argv[]) // the "index" command
return 1; return 1;
} }
if (prefix == 0) { if (prefix == 0) {
prefix = xmalloc(strlen(argv[optind]) + 4); prefix = malloc(strlen(argv[optind]) + 4);
strcpy(prefix, argv[optind]); strcpy(prefix, argv[optind]);
if (is_64) strcat(prefix, ".64"); if (is_64) strcat(prefix, ".64");
} }
str = (char*)xcalloc(strlen(prefix) + 10, 1); str = (char*)calloc(strlen(prefix) + 10, 1);
str2 = (char*)xcalloc(strlen(prefix) + 10, 1); str2 = (char*)calloc(strlen(prefix) + 10, 1);
str3 = (char*)xcalloc(strlen(prefix) + 10, 1); str3 = (char*)calloc(strlen(prefix) + 10, 1);
{ // nucleotide indexing { // nucleotide indexing
gzFile fp = xzopen(argv[optind], "r"); gzFile fp = xzopen(argv[optind], "r");

59
bwtsw2_aux.c
View File

@ -22,6 +22,11 @@ KSEQ_DECLARE(gzFile)
#define __left_lt(a, b) ((a).end > (b).end) #define __left_lt(a, b) ((a).end > (b).end)
KSORT_INIT(hit, bsw2hit_t, __left_lt) KSORT_INIT(hit, bsw2hit_t, __left_lt)
#ifdef USE_MALLOC_WRAPPERS
# include "malloc_wrap.h"
#endif
extern unsigned char nst_nt4_table[256]; extern unsigned char nst_nt4_table[256];
unsigned char nt_comp_table[256] = { unsigned char nt_comp_table[256] = {
@ -48,7 +53,7 @@ extern int bsw2_resolve_query_overlaps(bwtsw2_t *b, float mask_level);
bsw2opt_t *bsw2_init_opt() bsw2opt_t *bsw2_init_opt()
{ {
bsw2opt_t *o = (bsw2opt_t*)xcalloc(1, sizeof(bsw2opt_t)); bsw2opt_t *o = (bsw2opt_t*)calloc(1, sizeof(bsw2opt_t));
o->a = 1; o->b = 3; o->q = 5; o->r = 2; o->t = 30; o->a = 1; o->b = 3; o->q = 5; o->r = 2; o->t = 30;
o->bw = 50; o->bw = 50;
o->max_ins = 20000; o->max_ins = 20000;
@ -73,11 +78,11 @@ void bsw2_destroy(bwtsw2_t *b)
bwtsw2_t *bsw2_dup_no_cigar(const bwtsw2_t *b) bwtsw2_t *bsw2_dup_no_cigar(const bwtsw2_t *b)
{ {
bwtsw2_t *p; bwtsw2_t *p;
p = xcalloc(1, sizeof(bwtsw2_t)); p = calloc(1, sizeof(bwtsw2_t));
p->max = p->n = b->n; p->max = p->n = b->n;
if (b->n) { if (b->n) {
kroundup32(p->max); kroundup32(p->max);
p->hits = xcalloc(p->max, sizeof(bsw2hit_t)); p->hits = calloc(p->max, sizeof(bsw2hit_t));
memcpy(p->hits, b->hits, p->n * sizeof(bsw2hit_t)); memcpy(p->hits, b->hits, p->n * sizeof(bsw2hit_t));
} }
return p; return p;
@ -100,10 +105,10 @@ void bsw2_extend_left(const bsw2opt_t *opt, bwtsw2_t *b, uint8_t *_query, int lq
int8_t mat[25]; int8_t mat[25];
bwa_fill_scmat(opt->a, opt->b, mat); bwa_fill_scmat(opt->a, opt->b, mat);
query = xcalloc(lq, 1); query = calloc(lq, 1);
// sort according to the descending order of query end // sort according to the descending order of query end
ks_introsort(hit, b->n, b->hits); ks_introsort(hit, b->n, b->hits);
target = xcalloc(((lq + 1) / 2 * opt->a + opt->r) / opt->r + lq, 1); target = calloc(((lq + 1) / 2 * opt->a + opt->r) / opt->r + lq, 1);
// reverse _query // reverse _query
for (i = 0; i < lq; ++i) query[lq - i - 1] = _query[i]; for (i = 0; i < lq; ++i) query[lq - i - 1] = _query[i];
// core loop // core loop
@ -144,7 +149,7 @@ void bsw2_extend_rght(const bsw2opt_t *opt, bwtsw2_t *b, uint8_t *query, int lq,
int8_t mat[25]; int8_t mat[25];
bwa_fill_scmat(opt->a, opt->b, mat); bwa_fill_scmat(opt->a, opt->b, mat);
target = xcalloc(((lq + 1) / 2 * opt->a + opt->r) / opt->r + lq, 1); target = calloc(((lq + 1) / 2 * opt->a + opt->r) / opt->r + lq, 1);
for (i = 0; i < b->n; ++i) { for (i = 0; i < b->n; ++i) {
bsw2hit_t *p = b->hits + i; bsw2hit_t *p = b->hits + i;
int lt = ((lq - p->beg + 1) / 2 * opt->a + opt->r) / opt->r + lq; int lt = ((lq - p->beg + 1) / 2 * opt->a + opt->r) / opt->r + lq;
@ -192,7 +197,7 @@ static void gen_cigar(const bsw2opt_t *opt, int lq, uint8_t *seq[2], int64_t l_p
} }
#endif #endif
if (q->cigar && (beg != 0 || end < lq)) { // write soft clipping if (q->cigar && (beg != 0 || end < lq)) { // write soft clipping
q->cigar = xrealloc(q->cigar, 4 * (q->n_cigar + 2)); q->cigar = realloc(q->cigar, 4 * (q->n_cigar + 2));
if (beg != 0) { if (beg != 0) {
memmove(q->cigar + 1, q->cigar, q->n_cigar * 4); memmove(q->cigar + 1, q->cigar, q->n_cigar * 4);
q->cigar[0] = beg<<4 | 4; q->cigar[0] = beg<<4 | 4;
@ -223,7 +228,7 @@ static void merge_hits(bwtsw2_t *b[2], int l, int is_reverse)
int i; int i;
if (b[0]->n + b[1]->n > b[0]->max) { if (b[0]->n + b[1]->n > b[0]->max) {
b[0]->max = b[0]->n + b[1]->n; b[0]->max = b[0]->n + b[1]->n;
b[0]->hits = xrealloc(b[0]->hits, b[0]->max * sizeof(bsw2hit_t)); b[0]->hits = realloc(b[0]->hits, b[0]->max * sizeof(bsw2hit_t));
} }
for (i = 0; i < b[1]->n; ++i) { for (i = 0; i < b[1]->n; ++i) {
bsw2hit_t *p = b[0]->hits + b[0]->n + i; bsw2hit_t *p = b[0]->hits + b[0]->n + i;
@ -251,9 +256,9 @@ static bwtsw2_t *bsw2_aln1_core(const bsw2opt_t *opt, const bntseq_t *bns, uint8
_b = bsw2_core(bns, opt, query, target, pool); _b = bsw2_core(bns, opt, query, target, pool);
bwtl_destroy(query); bwtl_destroy(query);
for (k = 0; k < 2; ++k) { for (k = 0; k < 2; ++k) {
bb[k] = xcalloc(2, sizeof(void*)); bb[k] = calloc(2, sizeof(void*));
bb[k][0] = xcalloc(1, sizeof(bwtsw2_t)); bb[k][0] = calloc(1, sizeof(bwtsw2_t));
bb[k][1] = xcalloc(1, sizeof(bwtsw2_t)); bb[k][1] = calloc(1, sizeof(bwtsw2_t));
} }
for (k = 0; k < 2; ++k) { // separate _b into bb[2] based on the strand for (k = 0; k < 2; ++k) { // separate _b into bb[2] based on the strand
for (j = 0; j < _b[k]->n; ++j) { for (j = 0; j < _b[k]->n; ++j) {
@ -261,7 +266,7 @@ static bwtsw2_t *bsw2_aln1_core(const bsw2opt_t *opt, const bntseq_t *bns, uint8
p = bb[_b[k]->hits[j].is_rev][k]; p = bb[_b[k]->hits[j].is_rev][k];
if (p->n == p->max) { if (p->n == p->max) {
p->max = p->max? p->max<<1 : 8; p->max = p->max? p->max<<1 : 8;
p->hits = xrealloc(p->hits, p->max * sizeof(bsw2hit_t)); p->hits = realloc(p->hits, p->max * sizeof(bsw2hit_t));
} }
q = &p->hits[p->n++]; q = &p->hits[p->n++];
*q = _b[k]->hits[j]; *q = _b[k]->hits[j];
@ -340,7 +345,7 @@ static int fix_cigar(const bntseq_t *bns, bsw2hit_t *p, int n_cigar, uint32_t *c
uint32_t *cn; uint32_t *cn;
bwtint_t kk = 0; bwtint_t kk = 0;
nc = mq[0] = mq[1] = nlen[0] = nlen[1] = 0; nc = mq[0] = mq[1] = nlen[0] = nlen[1] = 0;
cn = xcalloc(n_cigar + 3, 4); cn = calloc(n_cigar + 3, 4);
x = coor; y = 0; x = coor; y = 0;
for (i = j = 0; i < n_cigar; ++i) { for (i = j = 0; i < n_cigar; ++i) {
int op = cigar[i]&0xf, ln = cigar[i]>>4; int op = cigar[i]&0xf, ln = cigar[i]>>4;
@ -398,9 +403,9 @@ static void write_aux(const bsw2opt_t *opt, const bntseq_t *bns, int qlen, uint8
if (b->n<<1 < b->max) { if (b->n<<1 < b->max) {
b->max = b->n; b->max = b->n;
kroundup32(b->max); kroundup32(b->max);
b->hits = xrealloc(b->hits, b->max * sizeof(bsw2hit_t)); b->hits = realloc(b->hits, b->max * sizeof(bsw2hit_t));
} }
b->aux = xcalloc(b->n, sizeof(bsw2aux_t)); b->aux = calloc(b->n, sizeof(bsw2aux_t));
// generate CIGAR // generate CIGAR
gen_cigar(opt, qlen, seq, bns->l_pac, pac, b, name); gen_cigar(opt, qlen, seq, bns->l_pac, pac, b, name);
// fix CIGAR, generate mapQ, and write chromosomal position // fix CIGAR, generate mapQ, and write chromosomal position
@ -559,7 +564,7 @@ static void bsw2_aln_core(bsw2seq_t *_seq, const bsw2opt_t *_opt, const bntseq_t
bsw2opt_t opt; bsw2opt_t opt;
bsw2global_t *pool = bsw2_global_init(); bsw2global_t *pool = bsw2_global_init();
bwtsw2_t **buf; bwtsw2_t **buf;
buf = xcalloc(_seq->n, sizeof(void*)); buf = calloc(_seq->n, sizeof(void*));
for (x = 0; x < _seq->n; ++x) { for (x = 0; x < _seq->n; ++x) {
bsw2seq1_t *p = _seq->seq + x; bsw2seq1_t *p = _seq->seq + x;
uint8_t *seq[2], *rseq[2]; uint8_t *seq[2], *rseq[2];
@ -570,10 +575,10 @@ static void bsw2_aln_core(bsw2seq_t *_seq, const bsw2opt_t *_opt, const bntseq_t
if (pool->max_l < l) { // then enlarge working space for aln_extend_core() if (pool->max_l < l) { // then enlarge working space for aln_extend_core()
int tmp = ((l + 1) / 2 * opt.a + opt.r) / opt.r + l; int tmp = ((l + 1) / 2 * opt.a + opt.r) / opt.r + l;
pool->max_l = l; pool->max_l = l;
pool->aln_mem = xrealloc(pool->aln_mem, (tmp + 2) * 24); pool->aln_mem = realloc(pool->aln_mem, (tmp + 2) * 24);
} }
// set seq[2] and rseq[2] // set seq[2] and rseq[2]
seq[0] = xcalloc(l * 4, 1); seq[0] = calloc(l * 4, 1);
seq[1] = seq[0] + l; seq[1] = seq[0] + l;
rseq[0] = seq[1] + l; rseq[1] = rseq[0] + l; rseq[0] = seq[1] + l; rseq[1] = rseq[0] + l;
// convert sequences to 2-bit representation // convert sequences to 2-bit representation
@ -586,7 +591,7 @@ static void bsw2_aln_core(bsw2seq_t *_seq, const bsw2opt_t *_opt, const bntseq_t
rseq[1][i] = c; rseq[1][i] = c;
} }
if (l - k < opt.t) { // too few unambiguous bases if (l - k < opt.t) { // too few unambiguous bases
buf[x] = xcalloc(1, sizeof(bwtsw2_t)); buf[x] = calloc(1, sizeof(bwtsw2_t));
free(seq[0]); continue; free(seq[0]); continue;
} }
// alignment // alignment
@ -618,7 +623,7 @@ static void bsw2_aln_core(bsw2seq_t *_seq, const bsw2opt_t *_opt, const bntseq_t
bsw2seq1_t *p = _seq->seq + x; bsw2seq1_t *p = _seq->seq + x;
uint8_t *seq[2]; uint8_t *seq[2];
int i; int i;
seq[0] = xmalloc(p->l * 2); seq[1] = seq[0] + p->l; seq[0] = malloc(p->l * 2); seq[1] = seq[0] + p->l;
for (i = 0; i < p->l; ++i) { for (i = 0; i < p->l; ++i) {
int c = nst_nt4_table[(int)p->seq[i]]; int c = nst_nt4_table[(int)p->seq[i]];
if (c >= 4) c = (int)(drand48() * 4); if (c >= 4) c = (int)(drand48() * 4);
@ -674,16 +679,16 @@ static void process_seqs(bsw2seq_t *_seq, const bsw2opt_t *opt, const bntseq_t *
int j; int j;
pthread_attr_init(&attr); pthread_attr_init(&attr);
pthread_attr_setdetachstate(&attr, PTHREAD_CREATE_JOINABLE); pthread_attr_setdetachstate(&attr, PTHREAD_CREATE_JOINABLE);
data = (thread_aux_t*)xcalloc(opt->n_threads, sizeof(thread_aux_t)); data = (thread_aux_t*)calloc(opt->n_threads, sizeof(thread_aux_t));
tid = (pthread_t*)xcalloc(opt->n_threads, sizeof(pthread_t)); tid = (pthread_t*)calloc(opt->n_threads, sizeof(pthread_t));
for (j = 0; j < opt->n_threads; ++j) { for (j = 0; j < opt->n_threads; ++j) {
thread_aux_t *p = data + j; thread_aux_t *p = data + j;
p->tid = j; p->_opt = opt; p->bns = bns; p->is_pe = is_pe; p->tid = j; p->_opt = opt; p->bns = bns; p->is_pe = is_pe;
p->pac = pac; p->target = target; p->pac = pac; p->target = target;
p->_seq = xcalloc(1, sizeof(bsw2seq_t)); p->_seq = calloc(1, sizeof(bsw2seq_t));
p->_seq->max = (_seq->n + opt->n_threads - 1) / opt->n_threads + 1; p->_seq->max = (_seq->n + opt->n_threads - 1) / opt->n_threads + 1;
p->_seq->n = 0; p->_seq->n = 0;
p->_seq->seq = xcalloc(p->_seq->max, sizeof(bsw2seq1_t)); p->_seq->seq = calloc(p->_seq->max, sizeof(bsw2seq1_t));
} }
for (i = 0; i < _seq->n; ++i) { // assign sequences to each thread for (i = 0; i < _seq->n; ++i) { // assign sequences to each thread
bsw2seq_t *p = data[(i>>is_pe)%opt->n_threads]._seq; bsw2seq_t *p = data[(i>>is_pe)%opt->n_threads]._seq;
@ -728,13 +733,13 @@ void bsw2_aln(const bsw2opt_t *opt, const bntseq_t *bns, bwt_t * const target, c
bsw2seq_t *_seq; bsw2seq_t *_seq;
bseq1_t *bseq; bseq1_t *bseq;
pac = xcalloc(bns->l_pac/4+1, 1); pac = calloc(bns->l_pac/4+1, 1);
for (l = 0; l < bns->n_seqs; ++l) for (l = 0; l < bns->n_seqs; ++l)
err_printf("@SQ\tSN:%s\tLN:%d\n", bns->anns[l].name, bns->anns[l].len); err_printf("@SQ\tSN:%s\tLN:%d\n", bns->anns[l].name, bns->anns[l].len);
err_fread_noeof(pac, 1, bns->l_pac/4+1, bns->fp_pac); err_fread_noeof(pac, 1, bns->l_pac/4+1, bns->fp_pac);
fp = xzopen(fn, "r"); fp = xzopen(fn, "r");
ks = kseq_init(fp); ks = kseq_init(fp);
_seq = xcalloc(1, sizeof(bsw2seq_t)); _seq = calloc(1, sizeof(bsw2seq_t));
if (fn2) { if (fn2) {
fp2 = xzopen(fn2, "r"); fp2 = xzopen(fn2, "r");
ks2 = kseq_init(fp2); ks2 = kseq_init(fp2);
@ -745,7 +750,7 @@ void bsw2_aln(const bsw2opt_t *opt, const bntseq_t *bns, bwt_t * const target, c
if (n > _seq->max) { if (n > _seq->max) {
_seq->max = n; _seq->max = n;
kroundup32(_seq->max); kroundup32(_seq->max);
_seq->seq = xrealloc(_seq->seq, _seq->max * sizeof(bsw2seq1_t)); _seq->seq = realloc(_seq->seq, _seq->max * sizeof(bsw2seq1_t));
} }
_seq->n = n; _seq->n = n;
for (i = 0; i < n; ++i) { for (i = 0; i < n; ++i) {

11
bwtsw2_chain.c
View File

@ -1,6 +1,9 @@
#include <stdio.h> #include <stdio.h>
#include "bwtsw2.h" #include "bwtsw2.h"
#include "utils.h"
#ifdef USE_MALLOC_WRAPPERS
# include "malloc_wrap.h"
#endif
typedef struct { typedef struct {
uint32_t tbeg, tend; uint32_t tbeg, tend;
@ -49,9 +52,9 @@ void bsw2_chain_filter(const bsw2opt_t *opt, int len, bwtsw2_t *b[2])
char *flag; char *flag;
// initialization // initialization
n[0] = b[0]->n; n[1] = b[1]->n; n[0] = b[0]->n; n[1] = b[1]->n;
z[0] = xcalloc(n[0] + n[1], sizeof(hsaip_t)); z[0] = calloc(n[0] + n[1], sizeof(hsaip_t));
z[1] = z[0] + n[0]; z[1] = z[0] + n[0];
chain[0] = xcalloc(n[0] + n[1], sizeof(hsaip_t)); chain[0] = calloc(n[0] + n[1], sizeof(hsaip_t));
for (k = j = 0; k < 2; ++k) { for (k = j = 0; k < 2; ++k) {
for (i = 0; i < b[k]->n; ++i) { for (i = 0; i < b[k]->n; ++i) {
bsw2hit_t *p = b[k]->hits + i; bsw2hit_t *p = b[k]->hits + i;
@ -74,7 +77,7 @@ void bsw2_chain_filter(const bsw2opt_t *opt, int len, bwtsw2_t *b[2])
} }
//for (k = 0; k < m[0]; ++k) printf("%d, [%d,%d), [%d,%d)\n", chain[0][k].chain, chain[0][k].tbeg, chain[0][k].tend, chain[0][k].qbeg, chain[0][k].qend); //for (k = 0; k < m[0]; ++k) printf("%d, [%d,%d), [%d,%d)\n", chain[0][k].chain, chain[0][k].tbeg, chain[0][k].tend, chain[0][k].qbeg, chain[0][k].qend);
// filtering // filtering
flag = xcalloc(m[0] + m[1], 1); flag = calloc(m[0] + m[1], 1);
ks_introsort(hsaip, m[0] + m[1], chain[0]); ks_introsort(hsaip, m[0] + m[1], chain[0]);
for (k = 1; k < m[0] + m[1]; ++k) { for (k = 1; k < m[0] + m[1]; ++k) {
hsaip_t *p = chain[0] + k; hsaip_t *p = chain[0] + k;

33
bwtsw2_core.c
View File

@ -7,7 +7,10 @@
#include "bwtsw2.h" #include "bwtsw2.h"
#include "bwt.h" #include "bwt.h"
#include "kvec.h" #include "kvec.h"
#include "utils.h"
#ifdef USE_MALLOC_WRAPPERS
# include "malloc_wrap.h"
#endif
typedef struct { typedef struct {
bwtint_t k, l; bwtint_t k, l;
@ -72,7 +75,7 @@ typedef struct __mempool_t {
inline static bsw2entry_p mp_alloc(mempool_t *mp) inline static bsw2entry_p mp_alloc(mempool_t *mp)
{ {
++mp->cnt; ++mp->cnt;
if (kv_size(mp->pool) == 0) return (bsw2entry_t*)xcalloc(1, sizeof(bsw2entry_t)); if (kv_size(mp->pool) == 0) return (bsw2entry_t*)calloc(1, sizeof(bsw2entry_t));
else return kv_pop(mp->pool); else return kv_pop(mp->pool);
} }
inline static void mp_free(mempool_t *mp, bsw2entry_p e) inline static void mp_free(mempool_t *mp, bsw2entry_p e)
@ -134,7 +137,7 @@ static void cut_tail(bsw2entry_t *u, int T, bsw2entry_t *aux)
if (u->n <= T) return; if (u->n <= T) return;
if (aux->max < u->n) { if (aux->max < u->n) {
aux->max = u->n; aux->max = u->n;
aux->array = (bsw2cell_t*)xrealloc(aux->array, aux->max * sizeof(bsw2cell_t)); aux->array = (bsw2cell_t*)realloc(aux->array, aux->max * sizeof(bsw2cell_t));
} }
a = (int*)aux->array; a = (int*)aux->array;
for (i = n = 0; i != u->n; ++i) for (i = n = 0; i != u->n; ++i)
@ -185,7 +188,7 @@ static void merge_entry(const bsw2opt_t * __restrict opt, bsw2entry_t *u, bsw2en
int i; int i;
if (u->n + v->n >= u->max) { if (u->n + v->n >= u->max) {
u->max = u->n + v->n; u->max = u->n + v->n;
u->array = (bsw2cell_t*)xrealloc(u->array, u->max * sizeof(bsw2cell_t)); u->array = (bsw2cell_t*)realloc(u->array, u->max * sizeof(bsw2cell_t));
} }
for (i = 0; i != v->n; ++i) { for (i = 0; i != v->n; ++i) {
bsw2cell_t *p = v->array + i; bsw2cell_t *p = v->array + i;
@ -203,7 +206,7 @@ static inline bsw2cell_t *push_array_p(bsw2entry_t *e)
{ {
if (e->n == e->max) { if (e->n == e->max) {
e->max = e->max? e->max<<1 : 256; e->max = e->max? e->max<<1 : 256;
e->array = (bsw2cell_t*)xrealloc(e->array, sizeof(bsw2cell_t) * e->max); e->array = (bsw2cell_t*)realloc(e->array, sizeof(bsw2cell_t) * e->max);
} }
return e->array + e->n; return e->array + e->n;
} }
@ -251,7 +254,7 @@ static void save_narrow_hits(const bwtl_t *bwtl, bsw2entry_t *u, bwtsw2_t *b1, i
if (p->G >= t && p->ql - p->qk + 1 <= IS) { // good narrow hit if (p->G >= t && p->ql - p->qk + 1 <= IS) { // good narrow hit
if (b1->max == b1->n) { if (b1->max == b1->n) {
b1->max = b1->max? b1->max<<1 : 4; b1->max = b1->max? b1->max<<1 : 4;
b1->hits = xrealloc(b1->hits, b1->max * sizeof(bsw2hit_t)); b1->hits = realloc(b1->hits, b1->max * sizeof(bsw2hit_t));
} }
q = &b1->hits[b1->n++]; q = &b1->hits[b1->n++];
q->k = p->qk; q->l = p->ql; q->k = p->qk; q->l = p->ql;
@ -280,7 +283,7 @@ int bsw2_resolve_duphits(const bntseq_t *bns, const bwt_t *bwt, bwtsw2_t *b, int
else if (p->G > 0) ++n; else if (p->G > 0) ++n;
} }
b->n = b->max = n; b->n = b->max = n;
b->hits = xcalloc(b->max, sizeof(bsw2hit_t)); b->hits = calloc(b->max, sizeof(bsw2hit_t));
for (i = j = 0; i < old_n; ++i) { for (i = j = 0; i < old_n; ++i) {
bsw2hit_t *p = old_hits + i; bsw2hit_t *p = old_hits + i;
if (p->l - p->k + 1 <= IS) { // the hit is no so repetitive if (p->l - p->k + 1 <= IS) { // the hit is no so repetitive
@ -400,9 +403,9 @@ bsw2global_t *bsw2_global_init()
{ {
bsw2global_t *pool; bsw2global_t *pool;
bsw2stack_t *stack; bsw2stack_t *stack;
pool = xcalloc(1, sizeof(bsw2global_t)); pool = calloc(1, sizeof(bsw2global_t));
stack = xcalloc(1, sizeof(bsw2stack_t)); stack = calloc(1, sizeof(bsw2stack_t));
stack->pool = (mempool_t*)xcalloc(1, sizeof(mempool_t)); stack->pool = (mempool_t*)calloc(1, sizeof(mempool_t));
pool->stack = (void*)stack; pool->stack = (void*)stack;
return pool; return pool;
} }
@ -462,13 +465,13 @@ bwtsw2_t **bsw2_core(const bntseq_t *bns, const bsw2opt_t *opt, const bwtl_t *ta
rhash = kh_init(qintv); rhash = kh_init(qintv);
init_bwtsw2(target, query, stack); init_bwtsw2(target, query, stack);
heap_size = opt->z; heap_size = opt->z;
heap = xcalloc(heap_size, sizeof(int)); heap = calloc(heap_size, sizeof(int));
// initialize the return struct // initialize the return struct
b = (bwtsw2_t*)xcalloc(1, sizeof(bwtsw2_t)); b = (bwtsw2_t*)calloc(1, sizeof(bwtsw2_t));
b->n = b->max = target->seq_len * 2; b->n = b->max = target->seq_len * 2;
b->hits = xcalloc(b->max, sizeof(bsw2hit_t)); b->hits = calloc(b->max, sizeof(bsw2hit_t));
b1 = (bwtsw2_t*)xcalloc(1, sizeof(bwtsw2_t)); b1 = (bwtsw2_t*)calloc(1, sizeof(bwtsw2_t));
b_ret = xcalloc(2, sizeof(void*)); b_ret = calloc(2, sizeof(void*));
b_ret[0] = b; b_ret[1] = b1; b_ret[0] = b; b_ret[1] = b1;
// initialize timer // initialize timer
getrusage(0, &last); getrusage(0, &last);

11
bwtsw2_pair.c
View File

@ -7,9 +7,12 @@
#include "bntseq.h" #include "bntseq.h"
#include "bwtsw2.h" #include "bwtsw2.h"
#include "kstring.h" #include "kstring.h"
#include "utils.h"
#include "ksw.h" #include "ksw.h"
#ifdef USE_MALLOC_WRAPPERS
# include "malloc_wrap.h"
#endif
#define MIN_RATIO 0.8 #define MIN_RATIO 0.8
#define OUTLIER_BOUND 2.0 #define OUTLIER_BOUND 2.0
#define MAX_STDDEV 4.0 #define MAX_STDDEV 4.0
@ -27,7 +30,7 @@ bsw2pestat_t bsw2_stat(int n, bwtsw2_t **buf, kstring_t *msg, int max_ins)
bsw2pestat_t r; bsw2pestat_t r;
memset(&r, 0, sizeof(bsw2pestat_t)); memset(&r, 0, sizeof(bsw2pestat_t));
isize = xcalloc(n, 8); isize = calloc(n, 8);
for (i = k = 0; i < n; i += 2) { for (i = k = 0; i < n; i += 2) {
bsw2hit_t *t[2]; bsw2hit_t *t[2];
int l; int l;
@ -116,7 +119,7 @@ void bsw2_pair1(const bsw2opt_t *opt, int64_t l_pac, const uint8_t *pac, const b
if (end > l_pac) end = l_pac; if (end > l_pac) end = l_pac;
if (end - beg < l_mseq) return; if (end - beg < l_mseq) return;
// generate the sequence // generate the sequence
seq = xmalloc(l_mseq + (end - beg)); seq = malloc(l_mseq + (end - beg));
ref = seq + l_mseq; ref = seq + l_mseq;
for (k = beg; k < end; ++k) for (k = beg; k < end; ++k)
ref[k - beg] = pac[k>>2] >> ((~k&3)<<1) & 0x3; ref[k - beg] = pac[k>>2] >> ((~k&3)<<1) & 0x3;
@ -195,7 +198,7 @@ void bsw2_pair(const bsw2opt_t *opt, int64_t l_pac, const uint8_t *pac, int n, b
a[which].flag |= BSW2_FLAG_RESCUED; a[which].flag |= BSW2_FLAG_RESCUED;
if (p[1]->max == 0) { if (p[1]->max == 0) {
p[1]->max = 1; p[1]->max = 1;
p[1]->hits = xmalloc(sizeof(bsw2hit_t)); p[1]->hits = malloc(sizeof(bsw2hit_t));
} }
p[1]->hits[0] = a[which]; p[1]->hits[0] = a[which];
p[1]->n = 1; p[1]->n = 1;

4
example.c
View File

@ -7,6 +7,10 @@
#include "kseq.h" // for the FASTA/Q parser #include "kseq.h" // for the FASTA/Q parser
KSEQ_DECLARE(gzFile) KSEQ_DECLARE(gzFile)
#ifdef USE_MALLOC_WRAPPERS
# include "malloc_wrap.h"
#endif
int main(int argc, char *argv[]) int main(int argc, char *argv[])
{ {
bwaidx_t *idx; bwaidx_t *idx;

6
is.c
View File

@ -27,6 +27,10 @@
#include <stdlib.h> #include <stdlib.h>
#include "utils.h" #include "utils.h"
#ifdef USE_MALLOC_WRAPPERS
# include "malloc_wrap.h"
#endif
typedef unsigned char ubyte_t; typedef unsigned char ubyte_t;
#define chr(i) (cs == sizeof(int) ? ((const int *)T)[i]:((const unsigned char *)T)[i]) #define chr(i) (cs == sizeof(int) ? ((const int *)T)[i]:((const unsigned char *)T)[i])
@ -205,7 +209,7 @@ int is_sa(const ubyte_t *T, int *SA, int n)
int is_bwt(ubyte_t *T, int n) int is_bwt(ubyte_t *T, int n)
{ {
int *SA, i, primary = 0; int *SA, i, primary = 0;
SA = (int*)xcalloc(n+1, sizeof(int)); SA = (int*)calloc(n+1, sizeof(int));
if (is_sa(T, SA, n)) err_fatal_simple("is_sa failed"); if (is_sa(T, SA, n)) err_fatal_simple("is_sa failed");

21
kbtree.h
View File

@ -31,7 +31,10 @@
#include <stdlib.h> #include <stdlib.h>
#include <string.h> #include <string.h>
#include <stdint.h> #include <stdint.h>
#include "utils.h"
#ifdef USE_MALLOC_WRAPPERS
# include "malloc_wrap.h"
#endif
typedef struct { typedef struct {
int32_t is_internal:1, n:31; int32_t is_internal:1, n:31;
@ -52,7 +55,7 @@ typedef struct {
kbtree_##name##_t *kb_init_##name(int size) \ kbtree_##name##_t *kb_init_##name(int size) \
{ \ { \
kbtree_##name##_t *b; \ kbtree_##name##_t *b; \
b = (kbtree_##name##_t*)xcalloc(1, sizeof(kbtree_##name##_t)); \ b = (kbtree_##name##_t*)calloc(1, sizeof(kbtree_##name##_t)); \
b->t = ((size - 4 - sizeof(void*)) / (sizeof(void*) + sizeof(key_t)) + 1) >> 1; \ b->t = ((size - 4 - sizeof(void*)) / (sizeof(void*) + sizeof(key_t)) + 1) >> 1; \
if (b->t < 2) { \ if (b->t < 2) { \
free(b); return 0; \ free(b); return 0; \
@ -61,7 +64,7 @@ typedef struct {
b->off_ptr = 4 + b->n * sizeof(key_t); \ b->off_ptr = 4 + b->n * sizeof(key_t); \
b->ilen = (4 + sizeof(void*) + b->n * (sizeof(void*) + sizeof(key_t)) + 3) >> 2 << 2; \ b->ilen = (4 + sizeof(void*) + b->n * (sizeof(void*) + sizeof(key_t)) + 3) >> 2 << 2; \
b->elen = (b->off_ptr + 3) >> 2 << 2; \ b->elen = (b->off_ptr + 3) >> 2 << 2; \
b->root = (kbnode_t*)xcalloc(1, b->ilen); \ b->root = (kbnode_t*)calloc(1, b->ilen); \
++b->n_nodes; \ ++b->n_nodes; \
return b; \ return b; \
} }
@ -70,7 +73,7 @@ typedef struct {
int i, max = 8; \ int i, max = 8; \
kbnode_t *x, **top, **stack = 0; \ kbnode_t *x, **top, **stack = 0; \
if (b) { \ if (b) { \
top = stack = (kbnode_t**)xcalloc(max, sizeof(kbnode_t*)); \ top = stack = (kbnode_t**)calloc(max, sizeof(kbnode_t*)); \
*top++ = (b)->root; \ *top++ = (b)->root; \
while (top != stack) { \ while (top != stack) { \
x = *--top; \ x = *--top; \
@ -79,7 +82,7 @@ typedef struct {
if (__KB_PTR(b, x)[i]) { \ if (__KB_PTR(b, x)[i]) { \
if (top - stack == max) { \ if (top - stack == max) { \
max <<= 1; \ max <<= 1; \
stack = (kbnode_t**)xrealloc(stack, max * sizeof(kbnode_t*)); \ stack = (kbnode_t**)realloc(stack, max * sizeof(kbnode_t*)); \
top = stack + (max>>1); \ top = stack + (max>>1); \
} \ } \
*top++ = __KB_PTR(b, x)[i]; \ *top++ = __KB_PTR(b, x)[i]; \
@ -173,7 +176,7 @@ typedef struct {
static void __kb_split_##name(kbtree_##name##_t *b, kbnode_t *x, int i, kbnode_t *y) \ static void __kb_split_##name(kbtree_##name##_t *b, kbnode_t *x, int i, kbnode_t *y) \
{ \ { \
kbnode_t *z; \ kbnode_t *z; \
z = (kbnode_t*)xcalloc(1, y->is_internal? b->ilen : b->elen); \ z = (kbnode_t*)calloc(1, y->is_internal? b->ilen : b->elen); \
++b->n_nodes; \ ++b->n_nodes; \
z->is_internal = y->is_internal; \ z->is_internal = y->is_internal; \
z->n = b->t - 1; \ z->n = b->t - 1; \
@ -211,7 +214,7 @@ typedef struct {
r = b->root; \ r = b->root; \
if (r->n == 2 * b->t - 1) { \ if (r->n == 2 * b->t - 1) { \
++b->n_nodes; \ ++b->n_nodes; \
s = (kbnode_t*)xcalloc(1, b->ilen); \ s = (kbnode_t*)calloc(1, b->ilen); \
b->root = s; s->is_internal = 1; s->n = 0; \ b->root = s; s->is_internal = 1; s->n = 0; \
__KB_PTR(b, s)[0] = r; \ __KB_PTR(b, s)[0] = r; \
__kb_split_##name(b, s, 0, r); \ __kb_split_##name(b, s, 0, r); \
@ -333,13 +336,13 @@ typedef struct {
#define __kb_traverse(key_t, b, __func) do { \ #define __kb_traverse(key_t, b, __func) do { \
int __kmax = 8; \ int __kmax = 8; \
__kbstack_t *__kstack, *__kp; \ __kbstack_t *__kstack, *__kp; \
__kp = __kstack = (__kbstack_t*)xcalloc(__kmax, sizeof(__kbstack_t)); \ __kp = __kstack = (__kbstack_t*)calloc(__kmax, sizeof(__kbstack_t)); \
__kp->x = (b)->root; __kp->i = 0; \ __kp->x = (b)->root; __kp->i = 0; \
for (;;) { \ for (;;) { \
while (__kp->x && __kp->i <= __kp->x->n) { \ while (__kp->x && __kp->i <= __kp->x->n) { \
if (__kp - __kstack == __kmax - 1) { \ if (__kp - __kstack == __kmax - 1) { \
__kmax <<= 1; \ __kmax <<= 1; \
__kstack = (__kbstack_t*)xrealloc(__kstack, __kmax * sizeof(__kbstack_t)); \ __kstack = (__kbstack_t*)realloc(__kstack, __kmax * sizeof(__kbstack_t)); \
__kp = __kstack + (__kmax>>1) - 1; \ __kp = __kstack + (__kmax>>1) - 1; \
} \ } \
(__kp+1)->i = 0; (__kp+1)->x = __kp->x->is_internal? __KB_PTR(b, __kp->x)[__kp->i] : 0; \ (__kp+1)->i = 0; (__kp+1)->x = __kp->x->is_internal? __KB_PTR(b, __kp->x)[__kp->i] : 0; \

11
khash.h
View File

@ -115,7 +115,10 @@ int main() {
#include <stdlib.h> #include <stdlib.h>
#include <string.h> #include <string.h>
#include <limits.h> #include <limits.h>
#include "utils.h"
#ifdef USE_MALLOC_WRAPPERS
# include "malloc_wrap.h"
#endif
/* compipler specific configuration */ /* compipler specific configuration */
@ -161,13 +164,13 @@ typedef khint_t khiter_t;
#endif #endif
#ifndef kcalloc #ifndef kcalloc
#define kcalloc(N,Z) xcalloc(N,Z) #define kcalloc(N,Z) calloc(N,Z)
#endif #endif
#ifndef kmalloc #ifndef kmalloc
#define kmalloc(Z) xmalloc(Z) #define kmalloc(Z) malloc(Z)
#endif #endif
#ifndef krealloc #ifndef krealloc
#define krealloc(P,Z) xrealloc(P,Z) #define krealloc(P,Z) realloc(P,Z)
#endif #endif
#ifndef kfree #ifndef kfree
#define kfree(P) free(P) #define kfree(P) free(P)

42
kopen.c
View File

@ -14,7 +14,9 @@
#include <sys/socket.h> #include <sys/socket.h>
#endif #endif
#include "utils.h" #ifdef USE_MALLOC_WRAPPERS
# include "malloc_wrap.h"
#endif
#ifdef _WIN32 #ifdef _WIN32
#define _KO_NO_NET #define _KO_NO_NET
@ -82,7 +84,7 @@ static int http_open(const char *fn)
// set ->http_host // set ->http_host
for (p = (char*)fn + 7; *p && *p != '/'; ++p); for (p = (char*)fn + 7; *p && *p != '/'; ++p);
l = p - fn - 7; l = p - fn - 7;
http_host = xcalloc(l + 1, 1); http_host = calloc(l + 1, 1);
strncpy(http_host, fn + 7, l); strncpy(http_host, fn + 7, l);
http_host[l] = 0; http_host[l] = 0;
for (q = http_host; *q && *q != ':'; ++q); for (q = http_host; *q && *q != ':'; ++q);
@ -91,21 +93,21 @@ static int http_open(const char *fn)
proxy = getenv("http_proxy"); proxy = getenv("http_proxy");
// set host, port and path // set host, port and path
if (proxy == 0) { if (proxy == 0) {
host = xstrdup(http_host); // when there is no proxy, server name is identical to http_host name. host = strdup(http_host); // when there is no proxy, server name is identical to http_host name.
port = xstrdup(*q? q : "80"); port = strdup(*q? q : "80");
path = xstrdup(*p? p : "/"); path = strdup(*p? p : "/");
} else { } else {
host = (strstr(proxy, "http://") == proxy)? xstrdup(proxy + 7) : xstrdup(proxy); host = (strstr(proxy, "http://") == proxy)? strdup(proxy + 7) : strdup(proxy);
for (q = host; *q && *q != ':'; ++q); for (q = host; *q && *q != ':'; ++q);
if (*q == ':') *q++ = 0; if (*q == ':') *q++ = 0;
port = xstrdup(*q? q : "80"); port = strdup(*q? q : "80");
path = xstrdup(fn); path = strdup(fn);
} }
/* connect; adapted from khttp_connect() in knetfile.c */ /* connect; adapted from khttp_connect() in knetfile.c */
l = 0; l = 0;
fd = socket_connect(host, port); fd = socket_connect(host, port);
buf = xcalloc(bufsz, 1); // FIXME: I am lazy... But in principle, 64KB should be large enough. buf = calloc(bufsz, 1); // FIXME: I am lazy... But in principle, 64KB should be large enough.
l += snprintf(buf + l, bufsz, "GET %s HTTP/1.0\r\nHost: %s\r\n\r\n", l += snprintf(buf + l, bufsz, "GET %s HTTP/1.0\r\nHost: %s\r\n\r\n",
path, http_host); path, http_host);
if (write_bytes(fd, buf, l) != 0) { if (write_bytes(fd, buf, l) != 0) {
@ -152,7 +154,7 @@ static int kftp_get_response(ftpaux_t *aux)
while (read(aux->ctrl_fd, &c, 1)) { // FIXME: this is *VERY BAD* for unbuffered I/O while (read(aux->ctrl_fd, &c, 1)) { // FIXME: this is *VERY BAD* for unbuffered I/O
if (n >= aux->max_response) { if (n >= aux->max_response) {
aux->max_response = aux->max_response? aux->max_response<<1 : 256; aux->max_response = aux->max_response? aux->max_response<<1 : 256;
aux->response = xrealloc(aux->response, aux->max_response); aux->response = realloc(aux->response, aux->max_response);
} }
aux->response[n++] = c; aux->response[n++] = c;
if (c == '\n') { if (c == '\n') {
@ -186,10 +188,10 @@ static int ftp_open(const char *fn)
for (p = (char*)fn + 6; *p && *p != '/'; ++p); for (p = (char*)fn + 6; *p && *p != '/'; ++p);
if (*p != '/') return 0; if (*p != '/') return 0;
l = p - fn - 6; l = p - fn - 6;
port = xstrdup("21"); port = strdup("21");
host = xcalloc(l + 1, 1); host = calloc(l + 1, 1);
strncpy(host, fn + 6, l); strncpy(host, fn + 6, l);
retr = xcalloc(strlen(p) + 8, 1); retr = calloc(strlen(p) + 8, 1);
sprintf(retr, "RETR %s\r\n", p); sprintf(retr, "RETR %s\r\n", p);
/* connect to ctrl */ /* connect to ctrl */
@ -241,8 +243,8 @@ static char **cmd2argv(const char *cmd)
for (i = beg + 1, argc = 0; i < end; ++i) for (i = beg + 1, argc = 0; i < end; ++i)
if (isspace(cmd[i]) && !isspace(cmd[i-1])) if (isspace(cmd[i]) && !isspace(cmd[i-1]))
++argc; ++argc;
argv = (char**)xcalloc(argc + 2, sizeof(void*)); argv = (char**)calloc(argc + 2, sizeof(void*));
argv[0] = str = (char*)xcalloc(end - beg + 1, 1); argv[0] = str = (char*)calloc(end - beg + 1, 1);
strncpy(argv[0], cmd + beg, end - beg); strncpy(argv[0], cmd + beg, end - beg);
for (i = argc = 1; i < end - beg; ++i) for (i = argc = 1; i < end - beg; ++i)
if (isspace(str[i])) str[i] = 0; if (isspace(str[i])) str[i] = 0;
@ -266,15 +268,15 @@ void *kopen(const char *fn, int *_fd)
koaux_t *aux = 0; koaux_t *aux = 0;
*_fd = -1; *_fd = -1;
if (strstr(fn, "http://") == fn) { if (strstr(fn, "http://") == fn) {
aux = xcalloc(1, sizeof(koaux_t)); aux = calloc(1, sizeof(koaux_t));
aux->type = KO_HTTP; aux->type = KO_HTTP;
aux->fd = http_open(fn); aux->fd = http_open(fn);
} else if (strstr(fn, "ftp://") == fn) { } else if (strstr(fn, "ftp://") == fn) {
aux = xcalloc(1, sizeof(koaux_t)); aux = calloc(1, sizeof(koaux_t));
aux->type = KO_FTP; aux->type = KO_FTP;
aux->fd = ftp_open(fn); aux->fd = ftp_open(fn);
} else if (strcmp(fn, "-") == 0) { } else if (strcmp(fn, "-") == 0) {
aux = xcalloc(1, sizeof(koaux_t)); aux = calloc(1, sizeof(koaux_t));
aux->type = KO_STDIN; aux->type = KO_STDIN;
aux->fd = STDIN_FILENO; aux->fd = STDIN_FILENO;
} else { } else {
@ -308,7 +310,7 @@ void *kopen(const char *fn, int *_fd)
exit(1); exit(1);
} else { /* parent process */ } else { /* parent process */
close(pfd[1]); close(pfd[1]);
aux = xcalloc(1, sizeof(koaux_t)); aux = calloc(1, sizeof(koaux_t));
aux->type = KO_PIPE; aux->type = KO_PIPE;
aux->fd = pfd[0]; aux->fd = pfd[0];
aux->pid = pid; aux->pid = pid;
@ -320,7 +322,7 @@ void *kopen(const char *fn, int *_fd)
*_fd = open(fn, O_RDONLY); *_fd = open(fn, O_RDONLY);
#endif #endif
if (*_fd >= 0) { if (*_fd >= 0) {
aux = xcalloc(1, sizeof(koaux_t)); aux = calloc(1, sizeof(koaux_t));
aux->type = KO_FILE; aux->type = KO_FILE;
aux->fd = *_fd; aux->fd = *_fd;
} }

21
kseq.h
View File

@ -31,7 +31,10 @@
#include <ctype.h> #include <ctype.h>
#include <string.h> #include <string.h>
#include <stdlib.h> #include <stdlib.h>
#include "utils.h"
#ifdef USE_MALLOC_WRAPPERS
# include "malloc_wrap.h"
#endif
#define KS_SEP_SPACE 0 // isspace(): \t, \n, \v, \f, \r #define KS_SEP_SPACE 0 // isspace(): \t, \n, \v, \f, \r
#define KS_SEP_TAB 1 // isspace() && !' ' #define KS_SEP_TAB 1 // isspace() && !' '
@ -51,9 +54,9 @@
#define __KS_BASIC(type_t, __bufsize) \ #define __KS_BASIC(type_t, __bufsize) \
static inline kstream_t *ks_init(type_t f) \ static inline kstream_t *ks_init(type_t f) \
{ \ { \
kstream_t *ks = (kstream_t*)xcalloc(1, sizeof(kstream_t)); \ kstream_t *ks = (kstream_t*)calloc(1, sizeof(kstream_t)); \
ks->f = f; \ ks->f = f; \
ks->buf = (unsigned char*)xmalloc(__bufsize); \ ks->buf = (unsigned char*)malloc(__bufsize); \
return ks; \ return ks; \
} \ } \
static inline void ks_destroy(kstream_t *ks) \ static inline void ks_destroy(kstream_t *ks) \
@ -121,7 +124,7 @@ typedef struct __kstring_t {
if (str->m - str->l < (size_t)(i - ks->begin + 1)) { \ if (str->m - str->l < (size_t)(i - ks->begin + 1)) { \
str->m = str->l + (i - ks->begin) + 1; \ str->m = str->l + (i - ks->begin) + 1; \
kroundup32(str->m); \ kroundup32(str->m); \
str->s = (char*)xrealloc(str->s, str->m); \ str->s = (char*)realloc(str->s, str->m); \
} \ } \
memcpy(str->s + str->l, ks->buf + ks->begin, i - ks->begin); \ memcpy(str->s + str->l, ks->buf + ks->begin, i - ks->begin); \
str->l = str->l + (i - ks->begin); \ str->l = str->l + (i - ks->begin); \
@ -133,7 +136,7 @@ typedef struct __kstring_t {
} \ } \
if (str->s == 0) { \ if (str->s == 0) { \
str->m = 1; \ str->m = 1; \
str->s = (char*)xcalloc(1, 1); \ str->s = (char*)calloc(1, 1); \
} else if (delimiter == KS_SEP_LINE && str->l > 1 && str->s[str->l-1] == '\r') --str->l; \ } else if (delimiter == KS_SEP_LINE && str->l > 1 && str->s[str->l-1] == '\r') --str->l; \
str->s[str->l] = '\0'; \ str->s[str->l] = '\0'; \
return str->l; \ return str->l; \
@ -152,7 +155,7 @@ typedef struct __kstring_t {
#define __KSEQ_BASIC(SCOPE, type_t) \ #define __KSEQ_BASIC(SCOPE, type_t) \
SCOPE kseq_t *kseq_init(type_t fd) \ SCOPE kseq_t *kseq_init(type_t fd) \
{ \ { \
kseq_t *s = (kseq_t*)xcalloc(1, sizeof(kseq_t)); \ kseq_t *s = (kseq_t*)calloc(1, sizeof(kseq_t)); \
s->f = ks_init(fd); \ s->f = ks_init(fd); \
return s; \ return s; \
} \ } \
@ -184,7 +187,7 @@ typedef struct __kstring_t {
if (c != '\n') ks_getuntil(ks, KS_SEP_LINE, &seq->comment, 0); /* read FASTA/Q comment */ \ if (c != '\n') ks_getuntil(ks, KS_SEP_LINE, &seq->comment, 0); /* read FASTA/Q comment */ \
if (seq->seq.s == 0) { /* we can do this in the loop below, but that is slower */ \ if (seq->seq.s == 0) { /* we can do this in the loop below, but that is slower */ \
seq->seq.m = 256; \ seq->seq.m = 256; \
seq->seq.s = (char*)xmalloc(seq->seq.m); \ seq->seq.s = (char*)malloc(seq->seq.m); \
} \ } \
while ((c = ks_getc(ks)) != -1 && c != '>' && c != '+' && c != '@') { \ while ((c = ks_getc(ks)) != -1 && c != '>' && c != '+' && c != '@') { \
if (c == '\n') continue; /* skip empty lines */ \ if (c == '\n') continue; /* skip empty lines */ \
@ -195,13 +198,13 @@ typedef struct __kstring_t {
if (seq->seq.l + 1 >= seq->seq.m) { /* seq->seq.s[seq->seq.l] below may be out of boundary */ \ if (seq->seq.l + 1 >= seq->seq.m) { /* seq->seq.s[seq->seq.l] below may be out of boundary */ \
seq->seq.m = seq->seq.l + 2; \ seq->seq.m = seq->seq.l + 2; \
kroundup32(seq->seq.m); /* rounded to the next closest 2^k */ \ kroundup32(seq->seq.m); /* rounded to the next closest 2^k */ \
seq->seq.s = (char*)xrealloc(seq->seq.s, seq->seq.m); \ seq->seq.s = (char*)realloc(seq->seq.s, seq->seq.m); \
} \ } \
seq->seq.s[seq->seq.l] = 0; /* null terminated string */ \ seq->seq.s[seq->seq.l] = 0; /* null terminated string */ \
if (c != '+') return seq->seq.l; /* FASTA */ \ if (c != '+') return seq->seq.l; /* FASTA */ \
if (seq->qual.m < seq->seq.m) { /* allocate memory for qual in case insufficient */ \ if (seq->qual.m < seq->seq.m) { /* allocate memory for qual in case insufficient */ \
seq->qual.m = seq->seq.m; \ seq->qual.m = seq->seq.m; \
seq->qual.s = (char*)xrealloc(seq->qual.s, seq->qual.m); \ seq->qual.s = (char*)realloc(seq->qual.s, seq->qual.m); \
} \ } \
while ((c = ks_getc(ks)) != -1 && c != '\n'); /* skip the rest of '+' line */ \ while ((c = ks_getc(ks)) != -1 && c != '\n'); /* skip the rest of '+' line */ \
if (c == -1) return -2; /* error: no quality string */ \ if (c == -1) return -2; /* error: no quality string */ \

9
ksort.h
View File

@ -57,7 +57,10 @@
#include <stdlib.h> #include <stdlib.h>
#include <string.h> #include <string.h>
#include "utils.h"
#ifdef USE_MALLOC_WRAPPERS
# include "malloc_wrap.h"
#endif
typedef struct { typedef struct {
void *left, *right; void *left, *right;
@ -73,7 +76,7 @@ typedef struct {
int curr, shift; \ int curr, shift; \
\ \
a2[0] = array; \ a2[0] = array; \
a2[1] = temp? temp : (type_t*)xmalloc(sizeof(type_t) * n); \ a2[1] = temp? temp : (type_t*)malloc(sizeof(type_t) * n); \
for (curr = 0, shift = 0; (1ul<<shift) < n; ++shift) { \ for (curr = 0, shift = 0; (1ul<<shift) < n; ++shift) { \
a = a2[curr]; b = a2[1-curr]; \ a = a2[curr]; b = a2[1-curr]; \
if (shift == 0) { \ if (shift == 0) { \
@ -183,7 +186,7 @@ typedef struct {
return; \ return; \
} \ } \
for (d = 2; 1ul<<d < n; ++d); \ for (d = 2; 1ul<<d < n; ++d); \
stack = (ks_isort_stack_t*)xmalloc(sizeof(ks_isort_stack_t) * ((sizeof(size_t)*d)+2)); \ stack = (ks_isort_stack_t*)malloc(sizeof(ks_isort_stack_t) * ((sizeof(size_t)*d)+2)); \
top = stack; s = a; t = a + (n-1); d <<= 1; \ top = stack; s = a; t = a + (n-1); d <<= 1; \
while (1) { \ while (1) { \
if (s < t) { \ if (s < t) { \

9
kstring.c
View File

@ -1,7 +1,10 @@
#include <stdarg.h> #include <stdarg.h>
#include <stdio.h> #include <stdio.h>
#include "kstring.h" #include "kstring.h"
#include "utils.h"
#ifdef USE_MALLOC_WRAPPERS
# include "malloc_wrap.h"
#endif
int ksprintf(kstring_t *s, const char *fmt, ...) int ksprintf(kstring_t *s, const char *fmt, ...)
{ {
@ -13,7 +16,7 @@ int ksprintf(kstring_t *s, const char *fmt, ...)
if (l + 1 > s->m - s->l) { if (l + 1 > s->m - s->l) {
s->m = s->l + l + 2; s->m = s->l + l + 2;
kroundup32(s->m); kroundup32(s->m);
s->s = (char*)xrealloc(s->s, s->m); s->s = (char*)realloc(s->s, s->m);
va_start(ap, fmt); va_start(ap, fmt);
l = vsnprintf(s->s + s->l, s->m - s->l, fmt, ap); l = vsnprintf(s->s + s->l, s->m - s->l, fmt, ap);
} }
@ -27,7 +30,7 @@ int ksprintf(kstring_t *s, const char *fmt, ...)
int main() int main()
{ {
kstring_t *s; kstring_t *s;
s = (kstring_t*)xcalloc(1, sizeof(kstring_t)); s = (kstring_t*)calloc(1, sizeof(kstring_t));
ksprintf(s, "abcdefg: %d", 100); ksprintf(s, "abcdefg: %d", 100);
printf("%s\n", s->s); printf("%s\n", s->s);
free(s); free(s);

17
kstring.h
View File

@ -3,7 +3,10 @@
#include <stdlib.h> #include <stdlib.h>
#include <string.h> #include <string.h>
#include "utils.h"
#ifdef USE_MALLOC_WRAPPERS
# include "malloc_wrap.h"
#endif
#ifndef kroundup32 #ifndef kroundup32
#define kroundup32(x) (--(x), (x)|=(x)>>1, (x)|=(x)>>2, (x)|=(x)>>4, (x)|=(x)>>8, (x)|=(x)>>16, ++(x)) #define kroundup32(x) (--(x), (x)|=(x)>>1, (x)|=(x)>>2, (x)|=(x)>>4, (x)|=(x)>>8, (x)|=(x)>>16, ++(x))
@ -22,7 +25,7 @@ static inline void ks_resize(kstring_t *s, size_t size)
if (s->m < size) { if (s->m < size) {
s->m = size; s->m = size;
kroundup32(s->m); kroundup32(s->m);
s->s = (char*)xrealloc(s->s, s->m); s->s = (char*)realloc(s->s, s->m);
} }
} }
@ -31,7 +34,7 @@ static inline int kputsn(const char *p, int l, kstring_t *s)
if (s->l + l + 1 >= s->m) { if (s->l + l + 1 >= s->m) {
s->m = s->l + l + 2; s->m = s->l + l + 2;
kroundup32(s->m); kroundup32(s->m);
s->s = (char*)xrealloc(s->s, s->m); s->s = (char*)realloc(s->s, s->m);
} }
memcpy(s->s + s->l, p, l); memcpy(s->s + s->l, p, l);
s->l += l; s->l += l;
@ -49,7 +52,7 @@ static inline int kputc(int c, kstring_t *s)
if (s->l + 1 >= s->m) { if (s->l + 1 >= s->m) {
s->m = s->l + 2; s->m = s->l + 2;
kroundup32(s->m); kroundup32(s->m);
s->s = (char*)xrealloc(s->s, s->m); s->s = (char*)realloc(s->s, s->m);
} }
s->s[s->l++] = c; s->s[s->l++] = c;
s->s[s->l] = 0; s->s[s->l] = 0;
@ -66,7 +69,7 @@ static inline int kputw(int c, kstring_t *s)
if (s->l + l + 1 >= s->m) { if (s->l + l + 1 >= s->m) {
s->m = s->l + l + 2; s->m = s->l + l + 2;
kroundup32(s->m); kroundup32(s->m);
s->s = (char*)xrealloc(s->s, s->m); s->s = (char*)realloc(s->s, s->m);
} }
for (x = l - 1; x >= 0; --x) s->s[s->l++] = buf[x]; for (x = l - 1; x >= 0; --x) s->s[s->l++] = buf[x];
s->s[s->l] = 0; s->s[s->l] = 0;
@ -83,7 +86,7 @@ static inline int kputuw(unsigned c, kstring_t *s)
if (s->l + l + 1 >= s->m) { if (s->l + l + 1 >= s->m) {
s->m = s->l + l + 2; s->m = s->l + l + 2;
kroundup32(s->m); kroundup32(s->m);
s->s = (char*)xrealloc(s->s, s->m); s->s = (char*)realloc(s->s, s->m);
} }
for (i = l - 1; i >= 0; --i) s->s[s->l++] = buf[i]; for (i = l - 1; i >= 0; --i) s->s[s->l++] = buf[i];
s->s[s->l] = 0; s->s[s->l] = 0;
@ -100,7 +103,7 @@ static inline int kputl(long c, kstring_t *s)
if (s->l + l + 1 >= s->m) { if (s->l + l + 1 >= s->m) {
s->m = s->l + l + 2; s->m = s->l + l + 2;
kroundup32(s->m); kroundup32(s->m);
s->s = (char*)xrealloc(s->s, s->m); s->s = (char*)realloc(s->s, s->m);
} }
for (x = l - 1; x >= 0; --x) s->s[s->l++] = buf[x]; for (x = l - 1; x >= 0; --x) s->s[s->l++] = buf[x];
s->s[s->l] = 0; s->s[s->l] = 0;

25
ksw.c
View File

@ -27,7 +27,10 @@
#include <stdint.h> #include <stdint.h>
#include <emmintrin.h> #include <emmintrin.h>
#include "ksw.h" #include "ksw.h"
#include "utils.h"
#ifdef USE_MALLOC_WRAPPERS
# include "malloc_wrap.h"
#endif
#ifdef __GNUC__ #ifdef __GNUC__
#define LIKELY(x) __builtin_expect((x),1) #define LIKELY(x) __builtin_expect((x),1)
@ -64,7 +67,7 @@ kswq_t *ksw_qinit(int size, int qlen, const uint8_t *query, int m, const int8_t
size = size > 1? 2 : 1; size = size > 1? 2 : 1;
p = 8 * (3 - size); // # values per __m128i p = 8 * (3 - size); // # values per __m128i
slen = (qlen + p - 1) / p; // segmented length slen = (qlen + p - 1) / p; // segmented length
q = (kswq_t*)xmalloc(sizeof(kswq_t) + 256 + 16 * slen * (m + 4)); // a single block of memory q = (kswq_t*)malloc(sizeof(kswq_t) + 256 + 16 * slen * (m + 4)); // a single block of memory
q->qp = (__m128i*)(((size_t)q + sizeof(kswq_t) + 15) >> 4 << 4); // align memory q->qp = (__m128i*)(((size_t)q + sizeof(kswq_t) + 15) >> 4 << 4); // align memory
q->H0 = q->qp + slen * m; q->H0 = q->qp + slen * m;
q->H1 = q->H0 + slen; q->H1 = q->H0 + slen;
@ -186,7 +189,7 @@ end_loop16:
if (n_b == 0 || (int32_t)b[n_b-1] + 1 != i) { // then append if (n_b == 0 || (int32_t)b[n_b-1] + 1 != i) { // then append
if (n_b == m_b) { if (n_b == m_b) {
m_b = m_b? m_b<<1 : 8; m_b = m_b? m_b<<1 : 8;
b = (uint64_t*)xrealloc(b, 8 * m_b); b = (uint64_t*)realloc(b, 8 * m_b);
} }
b[n_b++] = (uint64_t)imax<<32 | i; b[n_b++] = (uint64_t)imax<<32 | i;
} else if ((int)(b[n_b-1]>>32) < imax) b[n_b-1] = (uint64_t)imax<<32 | i; // modify the last } else if ((int)(b[n_b-1]>>32) < imax) b[n_b-1] = (uint64_t)imax<<32 | i; // modify the last
@ -289,7 +292,7 @@ end_loop8:
if (n_b == 0 || (int32_t)b[n_b-1] + 1 != i) { if (n_b == 0 || (int32_t)b[n_b-1] + 1 != i) {
if (n_b == m_b) { if (n_b == m_b) {
m_b = m_b? m_b<<1 : 8; m_b = m_b? m_b<<1 : 8;
b = (uint64_t*)xrealloc(b, 8 * m_b); b = (uint64_t*)realloc(b, 8 * m_b);
} }
b[n_b++] = (uint64_t)imax<<32 | i; b[n_b++] = (uint64_t)imax<<32 | i;
} else if ((int)(b[n_b-1]>>32) < imax) b[n_b-1] = (uint64_t)imax<<32 | i; // modify the last } else if ((int)(b[n_b-1]>>32) < imax) b[n_b-1] = (uint64_t)imax<<32 | i; // modify the last
@ -369,8 +372,8 @@ int ksw_extend(int qlen, const uint8_t *query, int tlen, const uint8_t *target,
int i, j, k, gapoe = gapo + gape, beg, end, max, max_i, max_j, max_gap, max_ie, gscore, max_off; int i, j, k, gapoe = gapo + gape, beg, end, max, max_i, max_j, max_gap, max_ie, gscore, max_off;
if (h0 < 0) h0 = 0; if (h0 < 0) h0 = 0;
// allocate memory // allocate memory
qp = xmalloc(qlen * m); qp = malloc(qlen * m);
eh = xcalloc(qlen + 1, 8); eh = calloc(qlen + 1, 8);
// generate the query profile // generate the query profile
for (k = i = 0; k < m; ++k) { for (k = i = 0; k < m; ++k) {
const int8_t *p = &mat[k * m]; const int8_t *p = &mat[k * m];
@ -461,7 +464,7 @@ static inline uint32_t *push_cigar(int *n_cigar, int *m_cigar, uint32_t *cigar,
if (*n_cigar == 0 || op != (cigar[(*n_cigar) - 1]&0xf)) { if (*n_cigar == 0 || op != (cigar[(*n_cigar) - 1]&0xf)) {
if (*n_cigar == *m_cigar) { if (*n_cigar == *m_cigar) {
*m_cigar = *m_cigar? (*m_cigar)<<1 : 4; *m_cigar = *m_cigar? (*m_cigar)<<1 : 4;
cigar = xrealloc(cigar, (*m_cigar) << 2); cigar = realloc(cigar, (*m_cigar) << 2);
} }
cigar[(*n_cigar)++] = len<<4 | op; cigar[(*n_cigar)++] = len<<4 | op;
} else cigar[(*n_cigar)-1] += len<<4; } else cigar[(*n_cigar)-1] += len<<4;
@ -477,9 +480,9 @@ int ksw_global(int qlen, const uint8_t *query, int tlen, const uint8_t *target,
if (n_cigar_) *n_cigar_ = 0; if (n_cigar_) *n_cigar_ = 0;
// allocate memory // allocate memory
n_col = qlen < 2*w+1? qlen : 2*w+1; // maximum #columns of the backtrack matrix n_col = qlen < 2*w+1? qlen : 2*w+1; // maximum #columns of the backtrack matrix
z = xmalloc(n_col * tlen); z = malloc(n_col * tlen);
qp = xmalloc(qlen * m); qp = malloc(qlen * m);
eh = xcalloc(qlen + 1, 8); eh = calloc(qlen + 1, 8);
// generate the query profile // generate the query profile
for (k = i = 0; k < m; ++k) { for (k = i = 0; k < m; ++k) {
const int8_t *p = &mat[k * m]; const int8_t *p = &mat[k * m];
@ -620,7 +623,7 @@ int main(int argc, char *argv[])
if (!forward_only) { // reverse if (!forward_only) { // reverse
if ((int)ksq->seq.m > max_rseq) { if ((int)ksq->seq.m > max_rseq) {
max_rseq = ksq->seq.m; max_rseq = ksq->seq.m;
rseq = (uint8_t*)xrealloc(rseq, max_rseq); rseq = (uint8_t*)realloc(rseq, max_rseq);
} }
for (i = 0, j = ksq->seq.l - 1; i < (int)ksq->seq.l; ++i, --j) for (i = 0, j = ksq->seq.l - 1; i < (int)ksq->seq.l; ++i, --j)
rseq[j] = ksq->seq.s[i] == 4? 4 : 3 - ksq->seq.s[i]; rseq[j] = ksq->seq.s[i] == 4? 4 : 3 - ksq->seq.s[i];

4
ksw.h
View File

@ -3,6 +3,10 @@
#include <stdint.h> #include <stdint.h>
#ifdef USE_MALLOC_WRAPPERS
# include "malloc_wrap.h"
#endif
#define KSW_XBYTE 0x10000 #define KSW_XBYTE 0x10000
#define KSW_XSTOP 0x20000 #define KSW_XSTOP 0x20000
#define KSW_XSUBO 0x40000 #define KSW_XSUBO 0x40000

12
kvec.h
View File

@ -50,6 +50,10 @@ int main() {
#include <stdlib.h> #include <stdlib.h>
#ifdef USE_MALLOC_WRAPPERS
# include "malloc_wrap.h"
#endif
#define kv_roundup32(x) (--(x), (x)|=(x)>>1, (x)|=(x)>>2, (x)|=(x)>>4, (x)|=(x)>>8, (x)|=(x)>>16, ++(x)) #define kv_roundup32(x) (--(x), (x)|=(x)>>1, (x)|=(x)>>2, (x)|=(x)>>4, (x)|=(x)>>8, (x)|=(x)>>16, ++(x))
#define kvec_t(type) struct { size_t n, m; type *a; } #define kvec_t(type) struct { size_t n, m; type *a; }
@ -60,7 +64,7 @@ int main() {
#define kv_size(v) ((v).n) #define kv_size(v) ((v).n)
#define kv_max(v) ((v).m) #define kv_max(v) ((v).m)
#define kv_resize(type, v, s) ((v).m = (s), (v).a = (type*)xrealloc((v).a, sizeof(type) * (v).m)) #define kv_resize(type, v, s) ((v).m = (s), (v).a = (type*)realloc((v).a, sizeof(type) * (v).m))
#define kv_copy(type, v1, v0) do { \ #define kv_copy(type, v1, v0) do { \
if ((v1).m < (v0).n) kv_resize(type, v1, (v0).n); \ if ((v1).m < (v0).n) kv_resize(type, v1, (v0).n); \
@ -71,19 +75,19 @@ int main() {
#define kv_push(type, v, x) do { \ #define kv_push(type, v, x) do { \
if ((v).n == (v).m) { \ if ((v).n == (v).m) { \
(v).m = (v).m? (v).m<<1 : 2; \ (v).m = (v).m? (v).m<<1 : 2; \
(v).a = (type*)xrealloc((v).a, sizeof(type) * (v).m); \ (v).a = (type*)realloc((v).a, sizeof(type) * (v).m); \
} \ } \
(v).a[(v).n++] = (x); \ (v).a[(v).n++] = (x); \
} while (0) } while (0)
#define kv_pushp(type, v) ((((v).n == (v).m)? \ #define kv_pushp(type, v) ((((v).n == (v).m)? \
((v).m = ((v).m? (v).m<<1 : 2), \ ((v).m = ((v).m? (v).m<<1 : 2), \
(v).a = (type*)xrealloc((v).a, sizeof(type) * (v).m), 0) \ (v).a = (type*)realloc((v).a, sizeof(type) * (v).m), 0) \
: 0), &(v).a[(v).n++]) : 0), &(v).a[(v).n++])
#define kv_a(type, v, i) (((v).m <= (size_t)(i)? \ #define kv_a(type, v, i) (((v).m <= (size_t)(i)? \
((v).m = (v).n = (i) + 1, kv_roundup32((v).m), \ ((v).m = (v).n = (i) + 1, kv_roundup32((v).m), \
(v).a = (type*)xrealloc((v).a, sizeof(type) * (v).m), 0) \ (v).a = (type*)realloc((v).a, sizeof(type) * (v).m), 0) \
: (v).n <= (size_t)(i)? (v).n = (i) + 1 \ : (v).n <= (size_t)(i)? (v).n = (i) + 1 \
: 0), (v).a[(i)]) : 0), (v).a[(i)])

57
malloc_wrap.c 100644
View File

@ -0,0 +1,57 @@
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <errno.h>
#ifdef USE_MALLOC_WRAPPERS
/* Don't wrap ourselves */
# undef USE_MALLOC_WRAPPERS
#endif
#include "malloc_wrap.h"
void *wrap_calloc(size_t nmemb, size_t size,
const char *file, unsigned int line, const char *func) {
void *p = calloc(nmemb, size);
if (NULL == p) {
fprintf(stderr,
"[%s] Failed to allocate %zd bytes at %s line %u: %s\n",
func, nmemb * size, file, line, strerror(errno));
exit(EXIT_FAILURE);
}
return p;
}
void *wrap_malloc(size_t size,
const char *file, unsigned int line, const char *func) {
void *p = malloc(size);
if (NULL == p) {
fprintf(stderr,
"[%s] Failed to allocate %zd bytes at %s line %u: %s\n",
func, size, file, line, strerror(errno));
exit(EXIT_FAILURE);
}
return p;
}
void *wrap_realloc(void *ptr, size_t size,
const char *file, unsigned int line, const char *func) {
void *p = realloc(ptr, size);
if (NULL == p) {
fprintf(stderr,
"[%s] Failed to allocate %zd bytes at %s line %u: %s\n",
func, size, file, line, strerror(errno));
exit(EXIT_FAILURE);
}
return p;
}
char *wrap_strdup(const char *s,
const char *file, unsigned int line, const char *func) {
char *p = strdup(s);
if (NULL == p) {
fprintf(stderr,
"[%s] Failed to allocate %zd bytes at %s line %u: %s\n",
func, strlen(s), file, line, strerror(errno));
exit(EXIT_FAILURE);
}
return p;
}

47
malloc_wrap.h 100644
View File

@ -0,0 +1,47 @@
#ifndef MALLOC_WRAP_H
#define MALLOC_WRAP_H
#include <stdlib.h> /* Avoid breaking the usual definitions */
#include <string.h>
#ifdef __cplusplus
extern "C" {
#endif
void *wrap_calloc(size_t nmemb, size_t size,
const char *file, unsigned int line, const char *func);
void *wrap_malloc(size_t size,
const char *file, unsigned int line, const char *func);
void *wrap_realloc(void *ptr, size_t size,
const char *file, unsigned int line, const char *func);
char *wrap_strdup(const char *s,
const char *file, unsigned int line, const char *func);
#ifdef __cplusplus
}
#endif
#ifdef USE_MALLOC_WRAPPERS
# ifdef calloc
# undef calloc
# endif
# define calloc(n, s) wrap_calloc( (n), (s), __FILE__, __LINE__, __func__)
# ifdef malloc
# undef malloc
# endif
# define malloc(s) wrap_malloc( (s), __FILE__, __LINE__, __func__)
# ifdef realloc
# undef realloc
# endif
# define realloc(p, s) wrap_realloc((p), (s), __FILE__, __LINE__, __func__)
# ifdef strdup
# undef strdup
# endif
# define strdup(s) wrap_strdup( (s), __FILE__, __LINE__, __func__)
#endif /* USE_MALLOC_WRAPPERS */
#endif /* MALLOC_WRAP_H */

18
pemerge.c
View File

@ -12,6 +12,10 @@
#include "utils.h" #include "utils.h"
KSEQ_DECLARE(gzFile) KSEQ_DECLARE(gzFile)
#ifdef USE_MALLOC_WRAPPERS
# include "malloc_wrap.h"
#endif
#define MAX_SCORE_RATIO 0.9f #define MAX_SCORE_RATIO 0.9f
#define MAX_ERR 8 #define MAX_ERR 8
@ -40,7 +44,7 @@ typedef struct {
pem_opt_t *pem_opt_init() pem_opt_t *pem_opt_init()
{ {
pem_opt_t *opt; pem_opt_t *opt;
opt = xcalloc(1, sizeof(pem_opt_t)); opt = calloc(1, sizeof(pem_opt_t));
opt->a = 5; opt->b = 4; opt->q = 2, opt->r = 17; opt->w = 20; opt->a = 5; opt->b = 4; opt->q = 2, opt->r = 17; opt->w = 20;
opt->T = opt->a * 10; opt->T = opt->a * 10;
opt->q_def = 20; opt->q_def = 20;
@ -58,8 +62,8 @@ int bwa_pemerge(const pem_opt_t *opt, bseq1_t x[2])
int i, xtra, l, l_seq, sum_q, ret = 0; int i, xtra, l, l_seq, sum_q, ret = 0;
kswr_t r; kswr_t r;
s[0] = xmalloc(x[0].l_seq); q[0] = xmalloc(x[0].l_seq); s[0] = malloc(x[0].l_seq); q[0] = malloc(x[0].l_seq);
s[1] = xmalloc(x[1].l_seq); q[1] = xmalloc(x[1].l_seq); s[1] = malloc(x[1].l_seq); q[1] = malloc(x[1].l_seq);
for (i = 0; i < x[0].l_seq; ++i) { for (i = 0; i < x[0].l_seq; ++i) {
int c = x[0].seq[i]; int c = x[0].seq[i];
s[0][i] = c < 0 || c > 127? 4 : c <= 4? c : nst_nt4_table[c]; s[0][i] = c < 0 || c > 127? 4 : c <= 4? c : nst_nt4_table[c];
@ -103,8 +107,8 @@ int bwa_pemerge(const pem_opt_t *opt, bseq1_t x[2])
l = x[0].l_seq - (r.tb - r.qb); // length to merge l = x[0].l_seq - (r.tb - r.qb); // length to merge
l_seq = x[0].l_seq + x[1].l_seq - l; l_seq = x[0].l_seq + x[1].l_seq - l;
seq = xmalloc(l_seq + 1); seq = malloc(l_seq + 1);
qual = xmalloc(l_seq + 1); qual = malloc(l_seq + 1);
memcpy(seq, s[0], x[0].l_seq); memcpy(seq + x[0].l_seq, &s[1][l], x[1].l_seq - l); memcpy(seq, s[0], x[0].l_seq); memcpy(seq + x[0].l_seq, &s[1][l], x[1].l_seq - l);
memcpy(qual, q[0], x[0].l_seq); memcpy(qual + x[0].l_seq, &q[1][l], x[1].l_seq - l); memcpy(qual, q[0], x[0].l_seq); memcpy(qual + x[0].l_seq, &q[1][l], x[1].l_seq - l);
for (i = 0, sum_q = 0; i < l; ++i) { for (i = 0, sum_q = 0; i < l; ++i) {
@ -174,7 +178,7 @@ static void process_seqs(const pem_opt_t *opt, int n_, bseq1_t *seqs, int64_t cn
int i, j, n = n_>>1<<1; int i, j, n = n_>>1<<1;
worker_t *w; worker_t *w;
w = xcalloc(opt->n_threads, sizeof(worker_t)); w = calloc(opt->n_threads, sizeof(worker_t));
for (i = 0; i < opt->n_threads; ++i) { for (i = 0; i < opt->n_threads; ++i) {
worker_t *p = &w[i]; worker_t *p = &w[i];
p->start = i; p->n = n; p->start = i; p->n = n;
@ -185,7 +189,7 @@ static void process_seqs(const pem_opt_t *opt, int n_, bseq1_t *seqs, int64_t cn
worker(w); worker(w);
} else { } else {
pthread_t *tid; pthread_t *tid;
tid = (pthread_t*)xcalloc(opt->n_threads, sizeof(pthread_t)); tid = (pthread_t*)calloc(opt->n_threads, sizeof(pthread_t));
for (i = 0; i < opt->n_threads; ++i) pthread_create(&tid[i], 0, worker, &w[i]); for (i = 0; i < opt->n_threads; ++i) pthread_create(&tid[i], 0, worker, &w[i]);
for (i = 0; i < opt->n_threads; ++i) pthread_join(tid[i], 0); for (i = 0; i < opt->n_threads; ++i) pthread_join(tid[i], 0);
free(tid); free(tid);

41
utils.c
View File

@ -264,47 +264,6 @@ int err_gzclose(gzFile file)
return ret; return ret;
} }
void *err_calloc(size_t nmemb, size_t size, const char *file, unsigned int line, const char *func)
{
void *p = calloc(nmemb, size);
if (NULL == p)
{
err_fatal(func, "Failed to allocate %zd bytes at %s line %u: %s\n", nmemb * size, file, line, strerror(errno));
}
return p;
}
void *err_malloc(size_t size, const char *file, unsigned int line, const char *func)
{
void *p = malloc(size);
if (NULL == p)
{
err_fatal(func, "Failed to allocate %zd bytes at %s line %u: %s\n", size, file, line, strerror(errno));
}
return p;
}
void *err_realloc(void *ptr, size_t size, const char *file, unsigned int line, const char *func)
{
void *p = realloc(ptr, size);
if (NULL == p)
{
err_fatal(func, "Failed to allocate %zd bytes at %s line %u: %s\n", size, file, line, strerror(errno));
}
return p;
}
char *err_strdup(const char *s, const char *file, unsigned int line, const char *func)
{
char *p = strdup(s);
if (NULL == p)
{
err_fatal(func, "Failed to allocate %zd bytes at %s line %u: %s\n", strlen(s), file, line, strerror(errno));
}
return p;
}
/********* /*********
* Timer * * Timer *
*********/ *********/

11
utils.h
View File

@ -48,12 +48,6 @@
#define xassert(cond, msg) if ((cond) == 0) _err_fatal_simple_core(__func__, msg) #define xassert(cond, msg) if ((cond) == 0) _err_fatal_simple_core(__func__, msg)
#define xcalloc(n, s) err_calloc( (n), (s), __FILE__, __LINE__, __func__)
#define xmalloc(s) err_malloc( (s), __FILE__, __LINE__, __func__)
#define xrealloc(p, s) err_realloc((p), (s), __FILE__, __LINE__, __func__)
#define xstrdup(s) err_strdup( (s), __FILE__, __LINE__, __func__)
typedef struct { typedef struct {
uint64_t x, y; uint64_t x, y;
} pair64_t; } pair64_t;
@ -91,11 +85,6 @@ extern "C" {
int err_fclose(FILE *stream); int err_fclose(FILE *stream);
int err_gzclose(gzFile file); int err_gzclose(gzFile file);
void *err_calloc(size_t nmemb, size_t size, const char *file, unsigned int line, const char *func);
void *err_malloc(size_t size, const char *file, unsigned int line, const char *func);
void *err_realloc(void *ptr, size_t size, const char *file, unsigned int line, const char *func);
char *err_strdup(const char *s, const char *file, unsigned int line, const char *func);
double cputime(); double cputime();
double realtime(); double realtime();