[Bug c++/16578] New: Seg fault @ Matrix multiplication algorithm. Recursion depth problem?

["schmitt at voidcast dot de" <gcc-bugzilla at gcc dot gnu dot org>]

Hi all, 
[SYSTEM: AMD Duron 700M,512M,Gentoo,gcc-3.3.3] 
[CODE main.cc:  
-------------- SNIPP ---------------- 
#include <stdio.h> 
#include <stdlib.h> 
#include <sys/time.h> 
 
 
template < typename atom, int n > struct StrassenMatrix { 
    atom data[n][n]; 
 
    StrassenMatrix < atom, n / 2 > operator() (int a, int b) { 
	StrassenMatrix < atom, n / 2 > result; 
	int offs_i = a * n / 2; 
	int offs_j = b * n / 2; 
	for (int i = 0; i < n / 2; i++) 
	    for (int j = 0; j < n / 2; j++) 
		 result.data[i][j] = data[i + offs_i][j + offs_j]; 
	 return result; 
    }  
     
    StrassenMatrix operator+(StrassenMatrix m) { 
	StrassenMatrix result; 
	for (int i = 0; i < n; i++) 
	    for (int j = 0; j < n; j++) 
		result.data[i][j] = data[i][j] + m.data[i][j]; 
	return result; 
    } 
 
    StrassenMatrix operator-(StrassenMatrix m) { 
	StrassenMatrix result; 
	for (int i = 0; i < n; i++) 
	    for (int j = 0; j < n; j++) 
		result.data[i][j] = data[i][j] - m.data[i][j]; 
	return result; 
    } 
 
    StrassenMatrix operator*(StrassenMatrix m2) { 
	StrassenMatrix result; 
	StrassenMatrix < atom, n / 2 > p1, p2, p3, p4, p5, p6, p7; 
	StrassenMatrix m1 = *this; 
	p1 = (m1(0, 1) - m1(1, 1)) * (m2(1, 0) + m2(1, 1)); 
	p2 = (m1(0, 0) + m1(1, 1)) * (m2(0, 0) + m2(1, 1)); 
	p3 = (m1(0, 0) - m1(1, 0)) * (m2(0, 0) + m2(0, 1)); 
	p4 = (m1(0, 0) + m1(0, 1)) * m2(1, 1); 
	p5 = m1(0, 0) * (m2(0, 1) - m2(1, 1)); 
	p6 = m1(1, 1) * (m2(1, 0) - m2(0, 0)); 
	p7 = (m1(1, 0) + m1(1, 1)) * m2(0, 0); 
	int nhalf = n / 2; 
	for (int i = 0; i < nhalf; i++) { 
	    for (int j = 0; j < nhalf; j++) { 
		result.data[i][j] =  
		    p1.data[i][j] + p2.data[i][j] - p4.data[i][j] +  
		    p6.data[i][j]; 
		result.data[i][j + nhalf] =  
		    p4.data[i][j] + p5.data[i][j]; 
		result.data[i + nhalf][j] =  
		    p6.data[i][j] + p7.data[i][j]; 
		result.data[i + nhalf][j + nhalf] = 
		    p2.data[i][j] - p3.data[i][j] + p5.data[i][j] - 
		    p7.data[i][j]; 
	    } 
	} 
	return result; 
    } 
 
    atom* elem(int i, int j) { 
	return &data[i][j]; 
    } 
}; 
 
template < typename atom > struct StrassenMatrix <atom, 1 > { 
    atom data[1][1]; 
 
    StrassenMatrix operator+(StrassenMatrix m) { 
	StrassenMatrix result; 
	 result.data[0][0] = data[0][0] + m.data[0][0]; 
	 return result; 
    }  
     
    StrassenMatrix operator-(StrassenMatrix m) { 
	StrassenMatrix result; 
	result.data[0][0] = data[0][0] - m.data[0][0]; 
	return result; 
    } 
 
    StrassenMatrix operator*(StrassenMatrix m) { 
	StrassenMatrix result; 
	result.data[0][0] = data[0][0] * m.data[0][0]; 
	return result; 
    } 
 
    atom* elem(int i, int j) { 
	return &data[i][j]; 
    } 
}; 
 
 
#define NUMTYPE double 
 
template<int dim> 
void test() { 
    // Initialisierung 
    StrassenMatrix < NUMTYPE, dim > strassen1, strassen2, strassenResult; 
    for (int i = 0; i < dim; i++) 
	for (int j = 0; j < dim; j++) { 
	    *strassen1.elem(i,j) = 
		(NUMTYPE) (10.0 * rand() / (RAND_MAX + 1.0)); 
	    *strassen2.elem(i,j) = 
		(NUMTYPE) (10.0 * rand() / (RAND_MAX + 1.0)); 
	} 
    // 5. Feld: Zeit fuer StrassenMatrix 
    strassenResult = strassen1 * strassen2; 
} 
 
int main(int argc, char **argv) 
{ 
    srand(0); 
    for (int i = 0; i < 5; i++) test<64>(); 
    for (int i = 0; i < 5; i++) test<128>(); 
    for (int i = 0; i < 5; i++) test<256>(); 
    for (int i = 0; i < 5; i++) test<512>(); 
    // ---*** HERE ***--- 
    // As soon as I uncomment the following line, g++ will seg fault.  
//    for (int i = 0; i < 5; i++) test<1024>(); 
    return 0; 
} 
-------------- SNIPP -------------- 
] 
 
======== 
SYMPTOM: 
When I do "g++ main.cc", and uncomment the line marked above with "HERE", 
g++ will seg fault. 
(The code above implements a matrix multiplication algorithm. In order to test 
its speed, I have it run with different matrix sizes. In order to see the 
magnitude of recursive creations of matrices, consider the --- not compiling 
-- 1024x1024 matrix multiplication: 
In the first step, the algorithm recurses to split up into four 512x512 
matrices each of which is then split up into four 256x256 matrices and so 
forth. So, the split-up-tree is 1:1024 - 4:512 - 16:256 - 64:128 - 256:64 - 
1024:32 - 4096:16 - 16384:8 - 65536:4 - 262144:2 - 1048576:1 
The --- compiling --- 512x512 matrix multiplication looks like this: 1:512 - 
4:256 - 16:128 - 64:64 - 256:32 - 1024:16 - 4096:8 - 16384:4 - 65536:2 - 
262144:1 
) 
 
REPRODUCEABILITY: 
When using gcc-3.3.3, this will always occur. 
Will not occur if compiled with gcc-2.95.3(or gcc-2.96? I am not totally 
sure; I have written and compiled this a while ago on a SuSE 8.2 installation, 
and it ran fine there.)  
Will not occur if compiled with gcc-3.3.1 (SuSE 9.0, just checked it) 
Will occur if compiled with gcc-3.4 (not totally sure if 3.4. Might have been 
3.5? tried this at a friend's machine. Was a Debian/unstable) 
 
COMMENT: 
I have no idea how a compiler works, I am just glad if it does what it should.  
So I won't speculate on the cause :-). 
For this problem though, I figure you don't need the output of  
"g++ -v -save-temps main.cc". If You do, I'd be glad to attach it. 
 
Greetings, 
Felix Schmitt.

-- 
           Summary: Seg fault @ Matrix multiplication algorithm. Recursion
                    depth problem?
           Product: gcc
           Version: 3.3.3
            Status: UNCONFIRMED
          Severity: minor
          Priority: P2
         Component: c++
        AssignedTo: unassigned at gcc dot gnu dot org
        ReportedBy: schmitt at voidcast dot de
                CC: gcc-bugs at gcc dot gnu dot org


http://gcc.gnu.org/bugzilla/show_bug.cgi?id=16578