/* Given a cutoff (suboptimal) score, this software computes all
   suboptimal pairs. */
/* The strength of this program is that it uses only space linear
   to the input (sequences' length) and output (all suboptimal pairs). */

/* Globally passed params and macros */

#include <stdio.h>
#include "sub.h"

/* CHECK_NODE - Check if (ii,jj) is a suboptimal grid point.
		If it is, create a new node and add it to the list. */
#define CHECK_NODE(ii,jj,cc,dd,ee,rr,ss,tt)				\
{   	if ((cc)+(rr)>=cut || (dd)+(ss)>=cut || (ee)+(tt)>=cut) {	\
	   /*Create a node */						\
	   np = (NODE *) ckalloc(sizeof(NODE));				\
	   np->i = (ii);						\
	   np->j = (jj);						\
	   np->c = (cc);						\
	   np->d = (dd);						\
	   np->e = (ee);						\
	   np->r = (rr);						\
	   np->s = (ss);						\
	   np->t = (tt);						\
	   np->link = NULL;						\
	   np->clink = NULL;						\
	   np->dlink = NULL;						\
	   np->elink = NULL;						\
	   np->state = 0;						\
									\
	   /* Add the node to the column list F */			\
	   if (F[(ii)]) {						\
		LastF[(ii)]->link = np;					\
		np->blink = LastF[(ii)];				\
	   } else {							\
		F[(ii)] = np;						\
		np->blink = NULL;					\
	   }								\
	   LastF[(ii)] = np;						\
	}								\
}

static long *CC, *DD;          /* Rightward score-only vectors */
static long *RR, *SS;          /* Leftward score-only vectors */
static long *ICC, *IDD;        /* Rightward score-only initial vectors */
static long *IRR, *ISS;        /* Leftward score-only initial vectors */
static long *IVC, *IVE;        /* Upward score-only initial vectors */
static long *IVR, *IVT;        /* Downward score-only initial vectors */
static long *C1, *C2;	       /* Temporary C vectors */
static long *D1, *D2;	       /* Temporary D vectors */
static long *E1, *E2;	       /* Temporary E vectors */
static long *R1, *R2;	       /* Temporary R vectors */
static long *S1, *S2;	       /* Temporary S vectors */
static long *T1, *T2;	       /* Temporary T vectors */

/* RIGHTWARD - forward computation */
#define RIGHTWARD					\
{ 	if ((c =   c   - m) > (e =   e   - h)) e = c; 	\
       	if ((c = CC[j] - m) > (d = DD[j] - h)) d = c; 	\
       	c = s + wa[BB[j]];				\
       	if (e > c) c = e;				\
       	if (d > c) c = d;				\
       	s = CC[j];					\
       	CC[j] = c;					\
       	DD[j] = d;					\
}

/* LEFTWARD - backward computation */
#define LEFTWARD			\
{ 	c = s + wa[BB[j+1]];		\
	d = SS[j]-h;			\
	if (e-m > c) c = e-m;		\
	if (d-g > c) c = d-g;		\
	if (c > d) d = c;		\
	if (c > (e = e-h)) e = c;	\
	s = RR[j];			\
	RR[j] = c;			\
	SS[j] = d;			\
}


/* sub_init - initializations for sub() */
void sub_init(int M, int N)
{
	register int i, j;
	long t;

	/* Compute ICC and IDD */
	ICC[0] = 0;
	t = -g;
	IDD[0] = MININT;
	for (j = 1; j <= N; ++j) {
		t = t-h;
		ICC[j] = t;
		IDD[j] = MININT;
	}
	/* Compute IVC and IVE */
	IVC[0] = 0;
	t = -g;
	IVE[0] = MININT;
	for (i = 1; i <= M; ++i) {
		t = t-h;
		IVC[i] = t;
		IVE[i] = MININT;
	}
	/* Compute IRR and ISS */
	IRR[N] = 0;
	ISS[N] = 0;
	t = -g;
	for (j = N-1; j >= 0; --j) {
		t = t-h;
		IRR[j] = t;
		ISS[j] = t;
	}
	/* Compute IVR and IVT */
	IVR[M] = 0;
	IVT[M] = 0;
	t = -g;
	for (i = M-1; i >= 0; --i) {
		t = t-h;
		IVR[i] = t;
		IVT[i] = t;
	}
}

/* sub - compute all suboptimal pairs in [I1,J1] x [I2,J2] */
/* 	 Initial vectors: BCC, BDD, BRR, BSS, BVC, BVE, BVR, BVT;
	 Lflag is TRUE if the leftmost column should be output in the
	 current subproblem; Bflag is TRUE if the lowest row should
	 be output in the current subproblem. */

int sub(int I1, int J1, int I2, int J2, long *BCC, long *BDD, long *BRR,
	 long *BSS, long *BVC, long *BVE, long *BVR, long *BVT, char Lflag,
	 char Bflag)
{
	int midI, midJ;		/* middle column and middle row */
	long *MCC, *MDD, *MRR, *MSS, *MVC, *MVE, *MVR, *MVT;
	long *MIS, *MJS, *LEFTS, *RIGHTS, *LOWS, *HIGHS;
	long s1, s2, s3, s4, s5, s6, s7, s8;
	int low, high, left, right; /* boundary nodes */
	char flag, flag1, flag2, flag3, flag4;
	int tI1, tJ1, tI2, tJ2;
	NODE *np;

	int i, j;
	long c, d, e, s, t, *wa;
	char lflag, bflag;
	char *ckalloc();

	/*
	printf("I1=%d, J1=%d, I2=%d, J2=%d, Lflag=%d, Bflag=%d\n",I1,J1,I2,J2,Lflag,Bflag);
	fflush(stdout);
	*/

	total_grid_pts += ((J2-J1+1) * (I2-I1+1));

	/* boundary cases: I1 = I2 or J1 = J2 */
	if (I1 >= I2) { /* one column */
	/*
	printf("I1=%d, J1=%d, I2=%d, J2=%d, Lflag=%d, Bflag=%d\n",I1,J1,I2,J2,Lflag,Bflag);
	*/
		if (Lflag) {
			/* temporary vectors */
			E1 = CC;
			T1 = RR;

			E1[J1] = BVE[I1];
			for (j=J1+1; j<=J2; ++j) {
				E1[j] = max(BCC[j-1]-m,E1[j-1]-h);
			}
			T1[J2] = BVT[I1];
			for (j=J2-1; j>=J1; --j) {
				T1[j] = max(BRR[j],T1[j+1]-h);
			}
			if (Bflag) {
			   CHECK_NODE(I1,J1,BCC[J1],BDD[J1],E1[J1],
					    BRR[J1],BSS[J1],T1[J1]);
			}
			for (j=J1+1; j<=J2; ++j) {
			   CHECK_NODE(I1,j,BCC[j],BDD[j],E1[j],
					    BRR[j],BSS[j],T1[j]);

			}
		}
		return;
	}
	if (J1 >= J2) { /* one row */
	/*
	printf("I1=%d, J1=%d, I2=%d, J2=%d, Lflag=%d, Bflag=%d\n",I1,J1,I2,J2,Lflag,Bflag);
	*/
		if (Bflag) {
			/* temporary vectors */
			j = (I2-I1+1)*sizeof(long);
			D1 = (long *) ckalloc(j) - I1;
			S1 = (long *) ckalloc(j) - I1;

			D1[I1] = BDD[J1];
			for (i=I1+1; i<=I2; ++i) {
				D1[i] = max(BVC[i-1]-m,D1[i-1]-h);
			}
			S1[I2] = BSS[J1];
			for (i=I2-1; i>=I1; --i) {
				S1[i] = max(BVR[i],S1[i+1]-h);
			}
			if (Lflag) {
			   CHECK_NODE(I1,J1,BVC[I1],D1[I1],BVE[I1],
					    BVR[I1],S1[I1],BVT[I1]);
			}
			for (i=I1+1; i<=I2; ++i) {
			   CHECK_NODE(i,J1,BVC[i],D1[i],BVE[i],
					    BVR[i],S1[i],BVT[i]);
			}
			free(D1+I1);
			free(S1+I1);
		}
		return;
	}

	/* boundary cases: I1+1 = I2 or J1+1 = J2 */
	if (I1+1 >= I2) { /* two columns */
	/*
	printf("I1=%d, J1=%d, I2=%d, J2=%d, Lflag=%d, Bflag=%d\n",I1,J1,I2,J2,Lflag,Bflag);
	*/
		/* Compute E1, C2, D2 and E2 */
		/* temporary vectors */
		E1 = CC;
		C2 = DD;
		D2 = RR;
		E2 = SS;

		E1[J1] = BVE[I1];
		C2[J1] = BVC[I2];
		D2[J1] = max(BCC[J1]-m, BDD[J1]-h);
		E2[J1] = BVE[I2];
		c = BVC[I2];
		e = BVE[I2];
      		wa = w[AA[I2]];
		for (j=J1+1; j<=J2; ++j) {
			E1[j] = max(BCC[j-1]-m, E1[j-1]-h);
 			if ((c =   c   - m) > (e =   e   - h)) e = c;
       			if ((c = BCC[j] - m) > (d = BDD[j] - h)) d = c;
       			c = BCC[j-1] + wa[BB[j]];
       			if (e > c) c = e;
       			if (d > c) c = d;
       			C2[j] = c;
       			D2[j] = d;
			E2[j] = e;
		}

		/* Compute T2, R1, S1 and T1 */
		/* temporary vectors */
		j = (J2-J1+1)*sizeof(long);
		R1 = (long *) ckalloc(j) - J1;
		S1 = (long *) ckalloc(j) - J1;
		T1 = (long *) ckalloc(j) - J1;
		T2 = (long *) ckalloc(j) - J1;

		T2[J2] = BVT[I2];
		R1[J2] = BVR[I1];
		S1[J2] = max(BVR[I1],BSS[J2]-h);
		T1[J2] = BVT[I1];

		wa = w[AA[I2]];
		e = T1[J2];
		for (j=J2-1; j>=J1; --j) {
			T2[j] = max(BRR[j],T2[j+1]-h);
 			c = BRR[j+1] + wa[BB[j+1]];
			d = BSS[j]-h;
			if (e-m > c) c = e-m;
			if (d-g > c) c = d-g;
			if (c > d) d = c;
			if (c > (e = e-h)) e = c;
			R1[j] = c;
			S1[j] = d;
			T1[j] = e;
		}

		if (Lflag && Bflag) {
			CHECK_NODE(I1,J1,BCC[J1],BDD[J1],E1[J1],
				   R1[J1],S1[J1],T1[J1]);
		}
		if (Lflag) {
			for (j=J1+1; j<=J2; ++j) {
			     CHECK_NODE(I1,j,BCC[j],BDD[j],E1[j],
				  R1[j],S1[j],T1[j]);
			}
		}
		if (Bflag) {
			CHECK_NODE(I2,J1,C2[J1],D2[J1],E2[J1],
				   BRR[J1],BSS[J1],T2[J1]);
		}
		for (j=J1+1; j<=J2; ++j) {
			CHECK_NODE(I2,j,C2[j],D2[j],E2[j],
				   BRR[j],BSS[j],T2[j]);
		}

		free(R1+J1);
		free(S1+J1);
		free(T1+J1);
		free(T2+J1);
		return;
	}
	if (J1+1 >= J2) { /* two rows */
	/*
	printf("I1=%d, J1=%d, I2=%d, J2=%d, Lflag=%d, Bflag=%d\n",I1,J1,I2,J2,Lflag,Bflag);
	*/
		/* temporary vectors */
		j = (I2-I1+1)*sizeof(long);
		D1 = (long *) ckalloc(j) - I1;
		R1 = (long *) ckalloc(j) - I1;
		S1 = (long *) ckalloc(j) - I1;
		T1 = (long *) ckalloc(j) - I1;
		C2 = (long *) ckalloc(j) - I1;
		D2 = (long *) ckalloc(j) - I1;
		E2 = (long *) ckalloc(j) - I1;
		S2 = (long *) ckalloc(j) - I1;

		/* Compute D1, C2, D2 and E2 */
		D1[I1] = BDD[J1];
		C2[I1] = BCC[J2];
		D2[I1] = BDD[J2];
		E2[I1] = max(BVC[I1]-m,BVE[I1]-h);
		c = C2[I1];
		d = D2[I1];
		for (i=I1+1; i<=I2; ++i) {
			D1[i] = max(BVC[i-1]-m,D1[i-1]-h);
       			if ((c = c - m) > (d = d - h)) d = c;
 			if ((c = BVC[i] - m) > (e = BVE[i] - h)) e = c;
       			c = BVC[i-1] + w[AA[i]][BB[J2]];
       			if (e > c) c = e;
       			if (d > c) c = d;
       			C2[i] = c;
       			D2[i] = d;
			E2[i] = e;
		}

		/* Compute S2, R1, S1 and T1 */
		S2[I2] = BSS[J2];
		R1[I2] = BRR[J1];
		S1[I2] = BSS[J1];
		T1[I2] = max(R1[I2],BVT[I2]-h);
		d = S1[I2];
		for (i=I2-1; i>=I1; --i) {
			S2[i] = max(BVR[i],S2[i+1]-h);
 			c = BVR[i+1] + w[AA[i+1]][BB[J2]];
			e = BVT[i]-h;
			if (d-m > c) c = d-m;
			if (e-g > c) c = e-g;
			if (c > e) e = c;
			if (c > (d = d-h)) d = c;
			R1[i] = c;
			S1[i] = d;
			T1[i] = e;
		}
		if (Lflag && Bflag) {
			CHECK_NODE(I1,J1,BVC[I1],D1[I1],BVE[I1],
				   R1[I1],S1[I1],T1[I1]);
		}
		if (Bflag) {
			for (i=I1+1; i<=I2; ++i) {
			    CHECK_NODE(i,J1,BVC[i],D1[i],BVE[i],
				   R1[i],S1[i],T1[i]);
			}
		}
		if (Lflag) {
			CHECK_NODE(I1,J2,C2[I1],D2[I1],E2[I1],
				   BVR[I1],S2[I1],BVT[I1]);
		}
		for (i=I1+1; i<=I2; ++i) {
			CHECK_NODE(i,J2,C2[i],D2[i],E2[i],
				   BVR[i],S2[i],BVT[i]);
		}

		free(D1+I1);
		free(R1+I1);
		free(S1+I1);
		free(T1+I1);
		free(C2+I1);
		free(D2+I1);
		free(E2+I1);
		free(S2+I1);
		return;
	}

	/* some initializations ...*/

	/* Initialize CC, DD, RR, SS */
	for (j = J1; j <= J2; ++j) {
		CC[j] = BCC[j];
		DD[j] = BDD[j];
		RR[j] = BRR[j];
		SS[j] = BSS[j];
	}

	/* Allocate space for middle column */
	j = (J2-J1+1) * sizeof(long);
	MCC = (long *) ckalloc(j) - J1;
	MDD = (long *) ckalloc(j) - J1;
	MRR = (long *) ckalloc(j) - J1;
	MSS = (long *) ckalloc(j) - J1;
	MIS = (long *) ckalloc(j) - J1;
	LEFTS = (long *) ckalloc(j) - J1;
	RIGHTS = (long *) ckalloc(j) - J1;

	/* Allocate space for middle row */
	j = (I2-I1+1) * sizeof(long);
	MVC = (long *) ckalloc(j) - I1;
	MVE = (long *) ckalloc(j) - I1;
	MVR = (long *) ckalloc(j) - I1;
	MVT = (long *) ckalloc(j) - I1;
	MJS = (long *) ckalloc(j) - I1;
	LOWS = (long *) ckalloc(j) - I1;
	HIGHS = (long *) ckalloc(j) - I1;

	/* Divide the problem by midI and midJ */

	midI = (I1+I2)/2;
	midJ = (J1+J2)/2;

	/* Rightward (forward) computation */

	/* Phase R1: compute [I1,J1] x [midI,J2] */
	/* Compute MVC[I1] and MVE[I1] */
	e = BVE[I1];
	for (j = J1+1; j <= midJ; ++j)
		e = max(BCC[j-1]-m, e-h);
	MVC[I1] = BCC[midJ];
	MVE[I1] = e;
	for (j = midJ+1; j <= J2; ++j)
		e = max(BCC[j-1]-m, e-h);
	HIGHS[I1] = max(BCC[J2]+BVR[I1],e+BVT[I1]);

	/* Compute from row I1+1 to row midI */
  	for (i = I1+1; i <= midI; ++i) {
     		s = CC[J1];
      		c = BVC[i];
      		CC[J1] = c;
		DD[J1] = max(s-m, DD[J1]-h);
      		e = BVE[i];
      		wa = w[AA[i]];
      		for (j = J1+1; j <= J2; ++j) {
			RIGHTWARD

			/* save c and e values on midJ */
			if (j == midJ) {
				MVC[i] = c;
				MVE[i] = e;
			}
		}
		HIGHS[i] = max(c+BVR[i],e+BVT[i]);
	}
	/* save c and d values on midI */
	for (j = J1; j <= J2; ++j) {
		MCC[j] = CC[j];
		MDD[j] = DD[j];
	}

	/* Phase R2: compute [midI+1,J1] x [I2,J2] */
  	for (i = midI+1; i <= I2; ++i) {
     		s = CC[J1];
      		c = BVC[i];
      		CC[J1] = c;
		DD[J1] = max(s-m, DD[J1]-h);
      		e = BVE[i];
      		wa = w[AA[i]];
      		for (j = J1+1; j <= J2; ++j) {
			RIGHTWARD

			if (j == midJ) {
				/* save c and e values on midJ */
				MVC[i] = c;
				MVE[i] = e;
			}
		}
		HIGHS[i] = max(c+BVR[i],e+BVT[i]);
	}

	/* Compute RIGHTS */
	for (j = J1; j <= J2; ++j)
		RIGHTS[j] = max(CC[j]+BRR[j],DD[j]+BSS[j]);

	/* Leftward (backward) computation */

	/* Phase L1: compute [midI,J1] x [I2,J2] */

	/* Compute MVR[I2] and MVT[I2] */
	e = BVT[I2];
	for (j = J2-1; j >= midJ; --j)
		e = max(BRR[j], e-h);
	MVR[I2] = BRR[midJ];
	MVT[I2] = e;
	for (j = midJ-1; j >= J1; --j)
		e = max(BRR[j], e-h);
	LOWS[I2] = max(BVC[I2]+BRR[J1], BVE[I2]+e);

	/* Compute from row I2-1 back to row midI */
	for (i = I2-1; i >= midI; --i) {
		s = RR[J2];
		c = BVR[i];
		RR[J2] = c;
		e = BVT[i];
		SS[J2] = max(c, SS[J2]-h);
		wa = w[AA[i+1]];
		for (j = J2-1; j >= J1; --j) {
			LEFTWARD

			/* save c and e values on midJ */
			if (j == midJ) {
				MVR[i] = c;
				MVT[i] = e;
			}
		}
		LOWS[i] = max(BVC[i]+c, BVE[i]+e);
	}

	/* save c and d values on midI */
	for (j = J1; j <= J2; ++j) {
		MRR[j] = RR[j];
		MSS[j] = SS[j];
	}

	/* Phase L2: compute [I1,J1] x [midI-1,J2] */
	for (i = midI-1; i >= I1; --i) {
		s = RR[J2];
		c = BVR[i];
		RR[J2] = c;
		e = BVT[i];
		SS[J2] = max(c, SS[J2]-h);
		wa = w[AA[i+1]];
		for (j = J2-1; j >= J1; --j) {
			LEFTWARD
			/* save c and e values on midJ */
			if (j == midJ) {
				MVR[i] = c;
				MVT[i] = e;
			}
		}
		LOWS[i] = max(BVC[i]+c, BVE[i]+e);
	}

	/* Compute LEFTS */
	for (j = J1; j <= J2; ++j)
		LEFTS[j] = max(BCC[j]+RR[j],BDD[j]+SS[j]);

	/* Compute MIS and MJS */
	for (j = J1; j <= J2; ++j)
		MIS[j] = max(MCC[j]+MRR[j],MDD[j]+MSS[j]);
	for (i = I1; i <= I2; ++i)
		MJS[i] = max(MVC[i]+MVR[i],MVE[i]+MVT[i]);

	/*
	for (j = J1; j <= J2; ++j)
		printf("MIS: j=%d, %d\n",j,MIS[j]);
	for (i = I1; i <= I2; ++i)
		printf("MJS: i=%d, %d\n",i,MJS[i]);
	*/

	/* Conquer: solve subproblems (up to 4) recursively */

	flag = FALSE;
	/* Subproblem 1: within [I1,J1] x [midI,midJ] */
	flag1 = FALSE;
	for (j = midJ; j >= J1; --j)
	   if (MIS[j] >= cut) {
		high = j;
		flag1 = TRUE;
		break;
	   }

	flag2 = FALSE;
	for (i = midI; i >= I1; --i)
	   if (MJS[i] >= cut) {
		right = i;
		flag2 = TRUE;
		break;
	   }
	/*
	printf("Subproblem 1: flag1 = %d, flag2 = %d\n",flag1,flag2);
	fflush(stdout);
	*/
	if (flag1 || flag2) {
		flag = TRUE;

		/* Save overlapped entries */
		s1 = BCC[midJ];
		s2 = BDD[midJ];
		s3 = MRR[midJ];
		s4 = MSS[midJ];
		s5 = BVC[midI];
		s6 = BVE[midI];
		s7 = MVR[midI];
		s8 = MVT[midI];

		tI1 = I1;
		tJ1 = J1;
		tI2 = midI;
		tJ2 = midJ;
		if (!flag1) {
		   tI2 = right;
		}
		if (!flag2) {
		   tJ2 = high;
		}
		flag3 = FALSE;
		for (j = J1; j <= midJ; ++j)
		   if (LEFTS[j] >= cut) {
			low = j;
			flag3 = TRUE;
			break;
		   }
		flag4 = FALSE;
		for (i = I1; i <= midI; ++i)
		   if (LOWS[i] >= cut) {
			left = i;
			flag4 = TRUE;
			break;
		   }
	/*
	printf("Subproblem 1: flag3 = %d, flag4 = %d\n",flag3,flag4);
	printf("high=%d, right=%d, low=%d, left=%d\n",high,right,low,left);
	fflush(stdout);
	*/
		lflag = Lflag;
		bflag = Bflag;
		if (!flag3) {
		   tI1 = left;
		}
		if (!flag4) {
		   tJ1 = low;
		}
		if (tI1 > I1) {
		   	lflag = TRUE;
			/* reset BCC and BDD */
			BCC[J1] = BVC[tI1];
			e = BVE[tI1];
			BDD[J1] = MININT;
			for (j = J1+1; j <= tJ2; ++j) {
				e = max(BCC[j-1]-m,e-h);
				BCC[j] = e;
				BDD[j] = MININT;
			}
		}
		if (tJ1 > J1) {
		   	bflag = TRUE;
			/* reset BVC and BVE */
			BVC[I1] = BCC[tJ1];
			d = BDD[tJ1];
			BVE[I1] = MININT;
			for (i = I1+1; i <= tI2; ++i) {
				d = max(BVC[i-1]-m,d-h);
				BVC[i] = d;
				BVE[i] = MININT;
			}
		}
		if (tI2 < midI) {
			/* reset MRR and MSS */
			MRR[midJ] = MVR[tI2];
			e = MVT[tI2];
			MSS[midJ] = MVR[tI2];
			for (j = midJ-1; j >= tJ1; --j) {
				MRR[j] = e-m;
				MSS[j] = MRR[j];
				e = e-h;
			}
		}
		if (tJ2 < midJ) {
			/* reset MVR and MVT */
			MVR[midI] = MRR[tJ2];
			d = MSS[tJ2];
			MVT[midI] = MRR[tJ2];
			for (i = midI-1; i >= tI1; --i) {
				MVR[i] = d-m;
				MVT[i] = MVR[i];
				d = d-h;
			}
		}
		sub(tI1,tJ1,tI2,tJ2,BCC,BDD,MRR,MSS,BVC,BVE,MVR,MVT,lflag,bflag);

		/* Restore overlapped entries */
		BCC[midJ] = s1;
		BDD[midJ] = s2;
		MRR[midJ] = s3;
		MSS[midJ] = s4;
		BVC[midI] = s5;
		BVE[midI] = s6;
		MVR[midI] = s7;
		MVT[midI] = s8;
	}

	/* Subproblem 2: within [I1,midJ] x [midI,J2] */
	flag1 = FALSE;
	for (j = J2; j >= midJ; --j)
	   if (MIS[j] >= cut) {
		high = j;
		flag1 = TRUE;
		break;
	   }

	flag2 = FALSE;
	for (i = midI; i >= I1; --i)
	   if (HIGHS[i] >= cut) {
		right = i;
		flag2 = TRUE;
		break;
	   }
	/*
	printf("Subproblem 2: flag1 = %d, flag2 = %d\n",flag1,flag2);
	fflush(stdout);
	*/
	if (flag1 || flag2) {
		flag = TRUE;

		/* Save overlapped entries */
		s1 = MVC[midI];
		s2 = MVE[midI];
		s3 = BVR[midI];
		s4 = BVT[midI];

		tI1 = I1;
		tJ1 = midJ;
		tI2 = midI;
		tJ2 = J2;
		if (!flag1) {
		   tI2 = right;
		}
		if (!flag2) {
		   tJ2 = high;
		}
		flag3 = FALSE;
		for (j = midJ; j <= J2; ++j)
		   if (LEFTS[j] >= cut) {
			low = j;
			flag3 = TRUE;
			break;
		   }
		flag4 = FALSE;
		for (i = I1; i <= midI; ++i)
		   if (MJS[i] >= cut) {
			left = i;
			flag4 = TRUE;
			break;
		   }
	/*
	printf("Subproblem 2: flag3 = %d, flag4 = %d\n",flag3,flag4);
	printf("high=%d, right=%d, low=%d, left=%d\n",high,right,low,left);
	fflush(stdout);
	*/
		lflag = Lflag;
		bflag = FALSE;
		if (!flag3) {
		   tI1 = left;
		}
		if (!flag4) {
		   tJ1 = low;
		}
		if (tI1 > I1) {
		   	lflag = TRUE;
			/* reset BCC and BDD */
			BCC[midJ] = MVC[tI1];
			e = MVE[tI1];
			BDD[midJ] = MININT;
			for (j = midJ+1; j <= tJ2; ++j) {
				e = max(BCC[j-1]-m,e-h);
				BCC[j] = e;
				BDD[j] = MININT;
			}
		}
		if (tJ1 > midJ) {
		   	bflag = TRUE;
			/* reset MVC and MVE */
			MVC[I1] = BCC[tJ1];
			d = BDD[tJ1];
			MVE[I1] = MININT;
			for (i = I1+1; i <= tI2; ++i) {
				d = max(MVC[i-1]-m,d-h);
				MVC[i] = d;
				MVE[i] = MININT;
			}
		}
		if (tI2 < midI) {
			/* reset MRR and MSS */
			MRR[J2] = BVR[tI2];
			e = BVT[tI2];
			MSS[J2] = BVR[tI2];
			for (j = J2-1; j >= tJ1; --j) {
				MRR[j] = e-m;
				MSS[j] = MRR[j];
				e = e-h;
			}
		}
		if (tJ2 < J2) {
			/* reset BVR and BVT */
			BVR[midI] = MRR[tJ2];
			d = MSS[tJ2];
			BVT[midI] = MRR[tJ2];
			for (i = midI-1; i >= tI1; --i) {
				BVR[i] = d-m;
				BVT[i] = BVR[i];
				d = d-h;
			}
		}
		sub(tI1,tJ1,tI2,tJ2,BCC,BDD,MRR,MSS,MVC,MVE,BVR,BVT,lflag,bflag);

		/* Restore overlapped entries */
		MVC[midI] = s1;
		MVE[midI] = s2;
		BVR[midI] = s3;
		BVT[midI] = s4;
	}


	/* Subproblem 3: within [midI,J1] x [I2,midJ] */
	flag1 = FALSE;
	for (j = midJ; j >= J1; --j)
	   if (RIGHTS[j] >= cut) {
		high = j;
		flag1 = TRUE;
		break;
	   }

	flag2 = FALSE;
	for (i = I2; i >= midI; --i)
	   if (MJS[i] >= cut) {
		right = i;
		flag2 = TRUE;
		break;
	   }
	/*
	printf("Subproblem 3: flag1 = %d, flag2 = %d\n",flag1,flag2);
	fflush(stdout);
	*/
	if (flag1 || flag2) {
		flag = TRUE;

		/* Save overlapped entries */
		s1 = MCC[midJ];
		s2 = MDD[midJ];
		s3 = BRR[midJ];
		s4 = BSS[midJ];

		tI1 = midI;
		tJ1 = J1;
		tI2 = I2;
		tJ2 = midJ;
		if (!flag1) {
		   tI2 = right;
		}
		if (!flag2) {
		   tJ2 = high;
		}
		flag3 = FALSE;
		for (j = J1; j <= midJ; ++j)
		   if (MIS[j] >= cut) {
			low = j;
			flag3 = TRUE;
			break;
		   }
		flag4 = FALSE;
		for (i = midI; i <= I2; ++i)
		   if (LOWS[i] >= cut) {
			left = i;
			flag4 = TRUE;
			break;
		   }
	/*
	printf("Subproblem 3: flag3 = %d, flag4 = %d\n",flag3,flag4);
	printf("high=%d, right=%d, low=%d, left=%d\n",high,right,low,left);
	fflush(stdout);
	*/
		lflag = FALSE;
		bflag = Bflag;
		if (!flag3) {
		   tI1 = left;
		}
		if (!flag4) {
		   tJ1 = low;
		}
		if (tI1 > midI) {
		   	lflag = TRUE;
			/* reset MCC and MDD */
			MCC[J1] = BVC[tI1];
			e = BVE[tI1];
			MDD[J1] = MININT;
			for (j = J1+1; j <= tJ2; ++j) {
				e = max(MCC[j-1]-m,e-h);
				MCC[j] = e;
				MDD[j] = MININT;
			}
		}
		if (tJ1 > J1) {
		   	bflag = TRUE;
			/* reset BVC and BVE */
			BVC[midI] = MCC[tJ1];
			d = MDD[tJ1];
			BVE[midI] = MININT;
			for (i = midI+1; i <= tI2; ++i) {
				d = max(BVC[i-1]-m,d-h);
				BVC[i] = d;
				BVE[i] = MININT;
			}
		}
		if (tI2 < I2) {
			/* reset BRR and BSS */
			BRR[midJ] = MVR[tI2];
			e = MVT[tI2];
			BSS[midJ] = MVR[tI2];
			for (j = midJ-1; j >= tJ1; --j) {
				BRR[j] = e-m;
				BSS[j] = BRR[j];
				e = e-h;
			}
		}
		if (tJ2 < midJ) {
			/* reset MVR and MVT */
			MVR[I2] = BRR[tJ2];
			d = BSS[tJ2];
			MVT[I2] = BRR[tJ2];
			for (i = I2-1; i >= tI1; --i) {
				MVR[i] = d-m;
				MVT[i] = MVR[i];
				d = d-h;
			}
		}
		sub(tI1,tJ1,tI2,tJ2,MCC,MDD,BRR,BSS,BVC,BVE,MVR,MVT,lflag,bflag);

		/* Restore overlapped entries */
		MCC[midJ] = s1;
		MDD[midJ] = s2;
		BRR[midJ] = s3;
		BSS[midJ] = s4;
	}

	/* Subproblem 4: within [midI,midJ] x [I2,J2] */
	flag1 = FALSE;
	for (j = J2; j >= midJ; --j)
	   if (RIGHTS[j] >= cut) {
		high = j;
		flag1 = TRUE;
		break;
	   }
	flag2 = FALSE;
	for (i = I2; i >= midI; --i)
	   if (HIGHS[i] >= cut) {
		right = i;
		flag2 = TRUE;
		break;
	   }
	/*
	printf("Subproblem 4: flag1 = %d, flag2 = %d\n",flag1,flag2);
	printf("midI=%d, midJ=%d\n",midI,midJ);
	fflush(stdout);
	*/
	if (flag1 || flag2) {
		flag = TRUE;

		tI1 = midI;
		tJ1 = midJ;
		tI2 = I2;
		tJ2 = J2;
		if (!flag1) {
		   tI2 = right;
		}
		if (!flag2) {
		   tJ2 = high;
		}
		flag3 = FALSE;
		for (j = midJ; j <= J2; ++j)
		   if (MIS[j] >= cut) {
			low = j;
			flag3 = TRUE;
			break;
		   }
		flag4 = FALSE;
		for (i = midI; i <= I2; ++i)
		   if (MJS[i] >= cut) {
			left = i;
			flag4 = TRUE;
			break;
		   }
	/*
	printf("Subproblem 4: flag3 = %d, flag4 = %d\n",flag3,flag4);
	printf("high=%d, right=%d, low=%d, left=%d\n",high,right,low,left);
	fflush(stdout);
	*/
		lflag = FALSE;
		bflag = FALSE;
		if (!flag3) {
		   tI1 = left;
		}
		if (!flag4) {
		   tJ1 = low;
		}
		if (tI1 > midI) {
		   	lflag = TRUE;
			/* reset MCC and MDD */
			MCC[midJ] = MVC[tI1];
			e = MVE[tI1];
			MDD[midJ] = MININT;
			/*
			printf("MCC:%d, MDD:%d, MVC:%d, MVE:%d\n",MCC[midJ],MDD[midJ],MVC[tI1],MVE[tI1]);
			printf("MVR:%d, MVT:%d\n",MVR[tI1],MVT[tI1]);
			*/
			for (j = midJ+1; j <= tJ2; ++j) {
				e = max(MCC[j-1]-m,e-h);
				MCC[j] = e;
				MDD[j] = MININT;
			}
		}
		if (tJ1 > midJ) {
		   	bflag = TRUE;
			/* reset MVC and MVE */
			MVC[midI] = MCC[tJ1];
			d = MDD[tJ1];
			MVE[midI] = MININT;
			for (i = midI+1; i <= tI2; ++i) {
				d = max(MVC[i-1]-m,d-h);
				MVC[i] = d;
				MVE[i] = MININT;
			}
		}
		if (tI2 < I2) {
			/* reset BRR and BSS */
			BRR[J2] = BVR[tI2];
			e = BVT[tI2];
			BSS[J2] = BVR[tI2];
			for (j = J2-1; j >= tJ1; --j) {
				BRR[j] = e-m;
				BSS[j] = BRR[j];
				e = e-h;
			}
		}
		if (tJ2 < J2) {
			/* reset BVR and BVT */
			BVR[I2] = BRR[tJ2];
			d = BSS[tJ2];
			BVT[I2] = BRR[tJ2];
			for (i = I2-1; i >= tI1; --i) {
				BVR[i] = d-m;
				BVT[i] = BVR[i];
				d = d-h;
			}
		}
		sub(tI1,tJ1,tI2,tJ2,MCC,MDD,BRR,BSS,MVC,MVE,BVR,BVT,lflag,bflag);
	}

	/* free middle column and row */
	free(MCC+J1);
	free(MDD+J1);
	free(MRR+J1);
	free(MSS+J1);
	free(MIS+J1);
	free(LEFTS+J1);
	free(RIGHTS+J1);
	free(MVC+I1);
	free(MVE+I1);
	free(MVR+I1);
	free(MVT+I1);
	free(MJS+I1);
	free(HIGHS+I1);
	free(LOWS+I1);
	if (flag) return 0;
	else return -1;
}

/*F_stat - print out some information about F */
F_stat(int M)
{
	int i, width, max_width_i, max_i, sp, total_sp, total_area;
	float avg_width, avg_sp;
	NODE *np;

	max_width_i = 0;
	total_sp = 0;
	total_area = 0;
	for (i=0; i<=M; ++i) {
		np = F[i];
		width = (LastF[i]->j - np->j + 1);
		if (width > max_width_i) {
			max_width_i = width;
			max_i = i;
		}
		total_area += width;
		sp = 0;
		while (np) {
			++sp;
			np = np->link;
		}
		total_sp += sp;
	}
	avg_width = (float) total_area / (float) (M+1);
	avg_sp = (float) total_sp / (float) (M+1);
	printf("max_i=%d, max_width_i=%d\n",max_i,max_width_i);
	printf("total_sp=%d, total_area=%d\n",total_sp,total_area);
	printf("avg_sp=%f, avg_width=%f\n",avg_sp,avg_width);
}

/* print_F - print all suboptimal gridpoints in their corresponding column
   	     in increasing row order. */
void print_F(int M)
{
	int i;
	NODE *np;

	for (i=0; i<=M; ++i) {
		printf(" ****** column %d ******\n",i);
		np = F[i];
		if (!np) printf("check column i=%d\n",i);
		while (np) {
		   printf("%d\n",np->j);
		   printf("row %d: %d %d %d %d %d %d\n",np->j,np->c,np->d,np->e,np->r,np->s,np->t);
		   np = np->link;
		}
	}
}

/* C_Construct - construct a column */
void C_Construct(int column)
{
	NODE  *p, *q;

	p = F[column];
	while (p) {
	   /* include edge e->c */
	   /*
	   if (p->e + p->r >= cut) {
		printf("e->c, %d\n",p->j);
	   }
	   */
	   /* include edge d->c */
	   /*
	   if (p->d + p->r >= cut) {
		printf("d->c, %d\n",p->j);
	   }
	   */
	   if (q = p->link) {
		if (p->j + 1 == q->j) {
		   if (p->e + q->t - h >= cut) {  /* include edge e->e */
			p->elink = q;
			/*
			printf("e->e, %d\n",p->j);
			*/
		   }
		   if (p->c + q->t - m >= cut) {  /* include edge c->e */
			p->elink = q;
			/*
			printf("c->e, %d\n",p->j);
			*/
		   }
		}
	   }
	   p = q;
	}
}

/* Column_link - link the current column to the previous column */
void Column_link(int column)
{
	NODE *p, *q;

	p = F[column];
	q = F[column-1];

	while (p) {
	   while (q && q->j < p->j - 1) q = q->link;
	   if (q && q->j == p->j - 1) { /* consider (column-1,q->j) -> (column,p->j) */
	      if (q->c + p->r + w[AA[column]][BB[p->j]] >= cut) {
		   /* include edge c->c (q->c to p->c) */
		   q->clink = p;
		   /*
		   printf("c->c, %d\n",p->j);
		   */
	      }
	      q = q->link;
	   }
	   if (q && q->j == p->j) {
	      if (q->d + p->s - h >= cut) { /* include edge d->d */
		   q->dlink = p;
		   /*
		   printf("d->d, %d\n",p->j);
		   */
	      }
	      if (q->c + p->s - m >= cut) { /* include edge c->d */
		   q->dlink = p;
		   /*
		   printf("c->d, %d\n",p->j);
		   */
	      }
	   }
	   p = p->link;
	}
}

/* FtoDAG - Transform F to a DAG */
void FtoDAG(int M)
{
	int i;

	/*
	printf("****** column %d ******\n",0);
	*/
	C_Construct(0);
	for (i=1; i<=M; ++i) {
		/*
		printf("****** column %d ******\n",i);
		*/
		C_Construct(i);
		Column_link(i);
	}
}

/* DAG_traverse - traverse the constructed suboptimal DAG */
void DAG_traverse(int M)
{
	int i;
	NODE *np;

	for (i=0; i<=M; ++i) {
		printf(" ****** column %d ******\n",i);
		np = F[i];
		if (!np) printf("check column i=%d\n",i);
		while (np) {
		   printf("%d",np->j);
		   if (np->clink) printf(" clink");
		   if (np->dlink) printf(" dlink");
		   if (np->elink) printf(" elink");
		   printf("\n");
		   if (!np->clink && !np->dlink && !np->elink)
			printf("(%d %d) is a dead end.\n",i,np->j);
		   np = np->link;
		}
	}
}

/* Interface */

int ALL_SUB(char *A, char *B, int M, int N, long W[][128], long G,
		long H, long CUT_OFF)
{
  int i, j, c;
  char *ckalloc();

  AA = A;			/* Setup global parameters */
  BB = B;
  w = W;
  g = G;
  h = H;
  m = g+h;
  cut = CUT_OFF;

  j = (N+1) * sizeof(long);
  CC = (long *) ckalloc(j);
  DD = (long *) ckalloc(j);
  RR = (long *) ckalloc(j);
  SS = (long *) ckalloc(j);
  ICC = (long *) ckalloc(j);
  IRR = (long *) ckalloc(j);
  IDD = (long *) ckalloc(j);
  ISS = (long *) ckalloc(j);
  j = (M+1) * sizeof(long);
  IVC = (long *) ckalloc(j);
  IVR = (long *) ckalloc(j);
  IVE = (long *) ckalloc(j);
  IVT = (long *) ckalloc(j);
  F = (NODE **) ckalloc((M+1)*sizeof(NODE *));
  for (i=0; i<=M; ++i) F[i] = NULL;
  LastF = (NODE **) ckalloc((M+1)*sizeof(NODE *));

  total_grid_pts = 0;
  sub_init(M,N);
  c = sub(0,0,M,N,ICC,IDD,IRR,ISS,IVC,IVE,IVR,IVT,TRUE,TRUE);

  free(CC); free(DD); free(RR); free(SS);
  free(ICC); free(IDD); free(IRR); free(ISS);
  free(IVC); free(IVR); free(IVE); free(IVT);

  printf("length of sequence 1: %d\n",M);
  printf("length of sequence 2: %d\n",N);
  printf("size = %d\n",M*N);
  printf("total grid points = %d\n",total_grid_pts);

  if (c == 0) {
  	/*
  	print_F(M);
  	*/
  	F_stat(M);
  	FtoDAG(M);
  	/*
  	DAG_traverse(M);
  	*/
	return 0;
  } else return -1;
}