Merge pull request ash211#8 from pgaurav/master

multiloop speedup changes

gtfold-mfe/include/energy.h

@@ -4,24 +4,20 @@
 #include "data.h"
 
 extern int *V; 
-extern int *VV; 
-extern int *VV1; 
 extern int *W; 
 extern int *VBI; 
 extern int *VM; 
 extern int **WM; 
-extern int *WMu; 
-extern int *WMl; 
+extern int **WMPrime; 
 extern int *indx; 
 
 
 #define V(i,j) V[indx[j]+i]
 #define VM(i,j) VM[indx[j]+i]
 #define WM(i,j) WM[i][j]
+#define WMPrime(i,j) WMPrime[i][j]
 #define WMU(i,j) WM[i][j]
 #define WML(i,j) WM[j][i]
-//#define WMU(i,j) WMu[indx[j]+i]
-//#define WML(i,j) WMl[indx[j]+i]
 #define VBI(i,j) VBI[indx[j]+i]
 
 #define auPen(i, j) ((( (i)==BASE_U || (j)==BASE_U ) && ( (i)==BASE_A || (i)==BASE_G || (j)==BASE_A || (j)==BASE_G )) ? auend : 0)

gtfold-mfe/src/algorithms.c

@@ -40,19 +40,19 @@ int calculate(int len, int nThreads) {
 	fprintf(stdout,"Thread count: %3d \n",omp_get_num_threads());
 #endif
 
-
 	for (b = TURN+1; b <= len-1; b++) {
 #ifdef _OPENMP
 #pragma omp parallel for private (i,j) schedule(guided)
 #endif
 		for (i = 1; i <= len - b; i++) {
 			j = i + b;
-			int flag = 0, newWM = INFINITY_; 
+			int newWM = INFINITY_; 
 			if (canPair(RNA[i], RNA[j])) {
-				flag = 1;
 				int eh = canHairpin(i,j)?eH(i,j):INFINITY_; //hair pin
-				int es = canStack(i,j)?eS(i,j)+getShapeEnergy(i)+getShapeEnergy(j)+V(i+1,j-1):INFINITY_; // stack
-				if (j-i > 6) { // Internal Loop BEGIN
+
+        int es = canStack(i,j)?eS(i,j)+getShapeEnergy(i)+getShapeEnergy(j)+V(i+1,j-1):INFINITY_; // stack
+
+        if (j-i > 6) { // Internal Loop BEGIN
 					int p=0, q=0;
 					int VBIij = INFINITY_;
 					for (p = i+1; p <= MIN(j-2-TURN,i+MAXLOOP+1) ; p++) {
@@ -67,18 +67,14 @@ int calculate(int len, int nThreads) {
 					}
 					VBI(i,j) = VBIij;
 					V(i,j) = V(i,j) + getShapeEnergy(i) + getShapeEnergy(j);
-
 				} // Internal Loop END
-				if (j-i > 10) { // Multi Loop BEGIN
-					int h;
-					int VMij, VMijd, VMidj, VMidjd;
-					VMij = VMijd = VMidj = VMidjd = INFINITY_;
-					for (h = i+TURN+1; h <= j-1-TURN; h++) { 
-						VMij = MIN(VMij, WMU(i+1,h-1) + WML(h,j-1)); 
-						VMidj = MIN(VMidj, WMU(i+2,h-1) + WML(h,j-1)); 
-						VMijd = MIN(VMijd, WMU(i+1,h-1) + WML(h,j-2)); 
-						VMidjd = MIN(VMidjd, WMU(i+2,h-1) + WML(h,j-2)); 
-					}
+
+        if (j-i > 10) { // Multi Loop BEGIN
+					int VMij =  WMPrime[i+1][j-1]; 
+					int VMidj = WMPrime[i+2][j-1]; 
+					int VMijd = WMPrime[i+1][j-2]; 
+					int VMidjd = WMPrime[i+2][j-2]; 
+
 					int d3 = canSS(j-1)?Ed3(i,j,j-1):INFINITY_;
 					int d5 = canSS(i+1)?Ed5(i,j,i+1):INFINITY_;
 					VMij = MIN(VMij, (VMidj + d5 +Ec)) ;
@@ -87,16 +83,21 @@ int calculate(int len, int nThreads) {
 					VMij = VMij + Ea + Eb + auPenalty(i,j);
 					VM(i,j) = canStack(i,j)?VMij:INFINITY_;
 				} // Multi Loop END
-				V(i,j) = MIN4(eh,es,VBI(i,j),VM(i,j));
+
+        V(i,j) = MIN4(eh,es,VBI(i,j),VM(i,j));
 			}
 			else V(i,j) = INFINITY_;
-			if (j-i > 4) { // WM BEGIN
+
+      if (j-i > 4) { // WM BEGIN
 				int h; 
 				for (h = i+TURN+1 ; h <= j-TURN-1; h++) {
 					//ZS: This sum corresponds to when i,j are NOT paired with each other.
 					//So we need to make sure only terms where i,j aren't pairing are considered. 
-					newWM = (!forcePair(i,j))?MIN(newWM, WMU(i,h-1) + WML(h,j)):newWM;
+					// Added auxillary storage WMPrime to speedup multiloop calculations
+          WMPrime[i][j] = MIN(WMPrime[i][j], WMU(i,h-1) + WML(h,j)); 
+					//newWM = (!forcePair(i,j))?MIN(newWM, WMU(i,h-1) + WML(h,j)):newWM;
 				}
+				newWM = (!forcePair(i,j))?MIN(newWM, WMPrime[i][j]):newWM;
 				newWM = MIN(V(i,j) + auPenalty(i,j) + Eb, newWM); 
 				newWM = canSS(i)?MIN(V(i+1,j) + Ed3(j,i+1,i) + auPenalty(i+1,j) + Eb + Ec, newWM):newWM; //i dangle
 				newWM = canSS(j)?MIN(V(i,j-1) + Ed5(j-1,i,j) + auPenalty(i,j-1) + Eb + Ec, newWM):newWM;  //j dangle
@@ -107,6 +108,7 @@ int calculate(int len, int nThreads) {
 			} // WM END
 		}
 	}
+
 	for (j = TURN+2; j <= len; j++) {
 		int i, Wj, Widjd, Wijd, Widj, Wij, Wim1;
 		Wj = INFINITY_;
@@ -123,4 +125,3 @@ int calculate(int len, int nThreads) {
 	}
 	return W[len];
 }
-

gtfold-mfe/src/energy.c

@@ -9,12 +9,11 @@
 #include "shapereader.h"
 
 int *V; 
-int *VV; 
-int *VV1; 
 int *W; 
 int *VBI; 
 int *VM; 
 int **WM; 
+int **WMPrime; 
 int *indx; 
 
 int alloc_flag = 0;
@@ -26,18 +25,6 @@ void create_tables(int len) {
 		exit(-1);
 	}
 
-	VV1 = (int *) malloc((len+ 1)*sizeof(int));
-	if (VV1 == NULL) {
-		perror("Cannot allocate variable 'V'");
-		exit(-1);
-	}
-
-	VV = (int *) malloc((len+ 1)*sizeof(int));
-	if (VV == NULL) {
-		perror("Cannot allocate variable 'V'");
-		exit(-1);
-	}
-
 	int i;
 	WM = (int **) malloc((len+1)* sizeof(int *));
 	if (WM == NULL) {
@@ -50,7 +37,20 @@ void create_tables(int len) {
 			perror("Cannot allocate variable 'WM[i]'");
 			exit(-1);
 		}   
+	}
+
+  WMPrime = (int **) malloc((len+1)* sizeof(int *));
+	if (WMPrime == NULL) {
+		perror("Cannot allocate variable 'WM'");
+		exit(-1);
 	}   
+	for (i = 0; i <= len; i++) {
+		WMPrime[i] = (int *)malloc((len+1)* sizeof(int));
+		if (WMPrime[i] == NULL) {
+			perror("Cannot allocate variable 'WM[i]'");
+			exit(-1);
+		}   
+	}
 
 	VM = (int *) malloc(((len+1)*len/2 + 1) * sizeof(int));
 	if (VM == NULL) {
@@ -87,15 +87,13 @@ void init_tables(int len) {
 
 	for (i = 0; i <= len; i++) {
 		W[i] = INFINITY_; 
-		VV[i] = INFINITY_;
-		VV1[i] = INFINITY_;
-		for (j = 0; j <= len; j++) 
+		for (j = 0; j <= len; j++) {
 			WM[i][j] = INFINITY_;
+			WMPrime[i][j] = INFINITY_;
+    }
 	}
 
-
 	LLL = (len)*(len+1)/2 + 1;
-
 	for (i = 0; i < LLL; i++) {
 		V[i] = INFINITY_;
 		VM[i] = INFINITY_;
@@ -116,11 +114,12 @@ void free_tables(int len) {
 		for (i = 0; i <= len; i++) free(WM[i]);
 		free(WM);
 
+		for (i = 0; i <= len; i++) free(WMPrime[i]);
+		free(WMPrime);
+
 		free(VM);
 		free(VBI);
 		free(V);
-		free(VV);
-		free(VV1);
 		free(W);
 	}
 }