迭代求子图同构

根据http://theory.stanford.edu/~virgi/cs267/lecture1.pdf论文第二部分的伪代码实现出来的C++实现代码，有改动，原文第5步没有还原。

#include<iostream>
#include<vector>
#include<time.h>
using namespace std;

#define MATRIX vector<vector<int>>
#define DEBUG

vector<int> getDegree(MATRIX m){
	int p = m.size();
	vector<int> degree(p,0);
	for (int i = 0; i < p; i++){
		for (int j = 0; j < p; j++){
			degree[i] += m[i][j];
		}
	}
	return degree;
}
MATRIX multi(MATRIX& M1, MATRIX& M2){
	int rowN = M1.size();
	int columnN = M2[0].size();
	MATRIX result;
	for (int i = 0; i < rowN; i++){
		result.push_back(vector<int>(columnN , 0));
		for (int j = 0; j < columnN; j++){
			for (int k = 0; k < M2.size(); k++){
				result[i][j] += M1[i][k] * M2[k][j];
			}
		}
	}
	return result;
}
MATRIX transposition(MATRIX& M){
	int rowN = M.size();
	int columnN = M[0].size();
	MATRIX result;
	for (int j = 0; j < columnN; j++){
		result.push_back(vector<int>(rowN));
		for (int i = 0; i < rowN; i++){
			result[j][i] = M[i][j];
		}
	}
	return result;
}
void printMatrix(MATRIX& M){
	for (MATRIX::iterator it = M.begin(); it != M.end(); it++){
		for (vector<int>::iterator itN = (*it).begin(); itN != (*it).end(); itN++){
			cout << *itN << " ";
		}
		cout << endl;
	}
}
bool checkM(MATRIX& M,MATRIX& MA,MATRIX& MB){
	MATRIX C = multi(M, transposition(multi(M, MB)));
	int rowN = MA.size();
	int columnN = MA[0].size();
	bool flag = true;
	for (int i = 0; i < rowN; i++){
		for (int j = 0; j < columnN; j++){
			if (MA[i][j] == 1 && C[i][j] != 1){
				flag = false;
				break;
			}
		}
		if (!flag)
			break;
	}
	return flag;
}

void ullmann(MATRIX &ma, MATRIX &mb){
	int pa = ma.size(), qa = ma[0].size();
	int pb = mb.size(), qb = mb[0].size();
	vector<int> degreeA = getDegree(ma);
	vector<int> degreeB = getDegree(mb);
	vector<vector<int>> M0;
	for (int i = 0; i < pa; i++){
		M0.push_back(vector<int>(pb));
		for (int j = 0; j < pb; j++){
			M0[i][j] = degreeB[j] >= degreeA[i] ? 1 : 0;	
		}
	}
	vector<int> F(pb, 0);
	vector<int> H(pa, -1);
	vector<MATRIX> M(pa);
	//step 1
	MATRIX CM = M0;
	int d = 0;
	H[0] = -1;
	int k = -1;
step2:
	bool existJ = false;
	for (int j = 0; j < pb; j++){
		if (CM[d][j] == 1 && F[j] == 0){
			existJ = true;
			break;
		}
	}
	if (existJ){
		M[d] = CM;
		k = d == 0 ? H[0] : -1;
	}
	else{
		goto step7;
	}
step3:
	k = k + 1;
	if (CM[d][k] == 0 || F[k] == 1)
		goto step3;
	H[d] = k;//
	for (int j = 0; j < pb; j++){
		if (j != k){
			CM[d][j] = 0;
		}
	}
step4:
	if (d < pa - 1 ){
		goto step6;
	}else{
		if (checkM(CM, ma, mb)){
#ifdef DEBUG
			//printMatrix(CM);
			static int sum = 0;//the number of answers
			cout << "---#" << ++sum << "---" << endl;
			for (int n = 0; n < H.size(); n++){
				cout << H[n]+1 << " ";
			}
			cout << endl;
#endif
		}
	}
	
step5:
	F[k] = 0;
	CM = M[d];//复原
	existJ = false;
	for (int j = k + 1; j < pb; j++){
		if (CM[d][j] == 1 && F[j] == 0){
			existJ = true;
			break;
		}
	}
	if (existJ){
		goto step3;
	}
	else{
		goto step7;
	}
step6:
	H[d] = k;
	F[k] = 1;
	//set F by H is right,or reset F[k] to 0 when search a new k in step 7.
	/*for (int i = 0; i < F.size(); i++)
		F[i] = 0;
	for (int i = 0; i <= d; i++){
		F[H[i]] = 1;
	}*/
	d = d + 1;
	goto step2;
step7:
	if (d == 0){
		goto END;
	}
	else{
		H[d] = -1;
		F[k] = 0;
		d = d - 1;
//		CM = M[d];
		k = H[d];
		goto step5;
	}
END:
#ifdef DEBUG
	cout <<"---end---"<< endl;
#endif
	return;
}

MATRIX getMatrix(int** a, int m, int n){
	MATRIX result;
	for (int i = 0; i < m; i++){
		result.push_back(vector<int>((int*)a + n*i, (int*)a + n*i + n));
	}
	return result;
}
int main(){
	int a[3][3] = { { 0, 1, 1 }, { 1, 0, 0 }, { 1, 0, 0 } };
	int b[4][4] = { { 0, 1, 1, 0 }, { 1, 0, 1, 1 }, { 1, 1, 0, 1 }, { 0, 1, 1, 0 } };
	MATRIX ma = getMatrix((int **)a, 3, 3);
	MATRIX mb = getMatrix((int **)b, 4, 4);
	clock_t t1 = clock();
	ullmann(ma, mb);
	clock_t t2 = clock();
	cout << "ull time is : " << t2 - t1 << endl;
	return 0;
}