由二维数组引发的一系列思考

关于C++中矩阵二维数组的一些思考

1.如何创建动态二维数组：（以int为例）

int** matrix = new int* [lines];
for(int i=0; i<lines; i++){
    matrix[i] = new int [rows];
}

由于内存地址的缘故，此时按行访问顺序会较快。

2.如何删除动态二维数组：

for(int i=0; i<lines; i++){
    delete[] matrix[i];
}
delete[] matrix;

如果动态删除二维数组的代码作为class的析构代码会发生什么呢？

~matrix(){
    //上述代码
}

让我们考虑以下情况：

//错误示范
class matrix(){//定义矩阵类matrix
public:
    int** m_matrix;
    int m_rows, m_lines;
public:
    matrix(int rows=2, int lines=2): m_rows(rows), m_lines(lines){
        m_matrix = new int* [m_lines];
        for(int i=0; i<lines; i++){
            m_matrix[i] = new int [m_rows];
        }
        for(int i=0; i<lines; i++)
            for(int j=0; j<rows; j++)
                m_matrix[i][j] = 0;
    }
    ~matrix(){
        for(int i=0; i<lines; i++){
            delete[] m_matrix[i];
        }
        delete[] m_matrix;
    }
}

matrix multiply(matrix A, matrix B){//矩阵乘法A*B
    matrix result(A.m_rows, B.m_lines);
    //……
    return result;
}
int main(){
    matrix A, B;
    //为A、B赋值，定义matrix C
    C = multiply(A,B);//此处会报错！
    return 0;
}

C=multiply(A,B)看似没什么问题，但是由于C++的默认拷贝构造是浅拷贝，只会进行简单的值复制。

而什么是值复制？值复制只会对类的数据成员进行复制，不包含类中指向类数据成员的指针！

所以会导致什么情况呢？

类对象C的m_matrix指针和multiply函数中result的m_matrix指针将会指向同一块内存地址。

multiply函数中的result已离开作用域，为了将result的结果赋值给C，默认拷贝构造函数将复制C的成员变量的值（忽略指针），然后调用拷贝构造函数对应的析构~matrix()，由于没复制矩阵matrix，导致接下来离开函数域调用的析构函数~matrix()要删除一个空的matrix，导致result被delete两次，由于只new了一次，此时程序将会异常终止。

过程为：constructor->copy constructor->destructor->destructor

是个很隐蔽的错误呢~

 

使用深拷贝以避免错误

深拷贝的实质是什么？就是复制指针。仍以上述矩阵为例：

class matrix{
    //……
    matrix(const matrix& M){
        m_rows = M.rows;
        m_lines = M.lines;//正式请使用getRows()和getLines()
        m_matrix = new int* [m_lines];//此处避免了指针冲突
        for(int i=0; i<m_lines; i++){
            m_matrix[i] = new int [m_rows];
        }
        for(int i=0; i<m_lines; i++){
            for(int j=0; j<m_rows; j++){
                m_matrix[i][j] = M.m_matrix[i][j];
            }    
        }
    }
}

可是现在问题又来了，如果写成

    matrix C(multiply(A, B));

编译器又会报错，提示左值不可以这么使用，左值又是怎么回事呢？

---

左值 lvalue 和 右值 rvalue

左值就是存储在内存中有确定地址的变量。右值定义为非左值，就是存储在内存中有不确定地址的变量。

我们想一想很容易就可以得出只有左值才能被赋值的结论，对右值这样一个寄存器中的临时存储赋值是没有意义的，并且也不被允许。

在这样的前提下，我们来讨论以下几种情况：

1.能对函数返回值赋值吗

int x = 10;
int foo1(){
    return x;
}
int& foo2(){
    return x;
}
int main(){
    int v = 2;
    foo1(v) = 0;//错误，foo1(v)返回一个右值
    foo2(v) = 0;//正确，foo2(v)返回一个左值
    return 0;
}

代码说明当返回值为左值时可以被赋值，右值不可

2.什么样的左值可以被赋值

由于const关键字的存在，类似const int v=0; v=1;这种代码肯定是错的。

所以可修改的左值才可以被赋值。

3.左右值怎么转换

首先，当我们将左值放在等号右边时，左值将隐式地转换为右值，就是从内存复制到寄存器缓存中去了，例如：

int a = 0, b = 1, c = 2;
int c = a + b;//(a+b)就是右值

显式的包括引用（&）等，例如：

int value;
int* p = &value;//正确
int* q = &(value + 1);//错误，不能对右值取地址

右值转为左值可以用解引用（*）等，例如：

int value[2] = {0, 1};
int* p = &value[0];
*(p + 1) = 5;//右值p+1转换为左值

4.关于引用&

C++中&除了作为取地址符，还作为引用类型标识符，注意int&是一个整体，不能分开再合上看含义

int& b = a; 等价于int* const b = & a;，const表明b始终指向a，b改变将引起a改变。

注意下面写法中存在的问题：

int cun=0;
int* p=&cun;
int& mura=*p;//等价于int* const mura=*p（解引用）=cun
int bun=1;
p=&bun;//mura是const变量而p却可以改变，这与引用的初衷相违背

此时mura还是引用cun不变，因为引用只能也必须初始化，而不能之后再去赋值。

那么引用作为函数变量时会发生什么？

简单的说就是将真实的变量传进了函数，而不是复制出的临时变量，所以函数中参变量改变将会引起传入变量的改变。

如果不希望函数对引用变量作出修改，建议写成：

void add2(const int& a, const int& b){……}

如果形参与实际传入的参数类型不一致（例如long 和 int），仅在形参声明为const的情况下，会生成临时变量代替实参，也就是实参不会被改变，未声明const不一致直接报错。

扩展：右值引用&&

右值不可取地址，但是我们可以令一个变量关联右值，那么我们就可以通过这个变量来获取右值了。

long foo(){return 12;}
int x=1, y=2;
int&& r=x+y;
double&& s=3.1415926;
long&& t=foo();

取用时，写r，s，t即可取值，也可&r取用r的地址

右值引用在一些参变量占很大内存时用处极大，因为在跳出函数体时，返回值将被复制一份用于赋值，然后再去销毁，浪费极大空间时间。

---

回到刚才的代码

class matrix{
    //……
    matrix(const matrix& M){……}
}
……
    matrix C(multiply(A, B));

显然犯了将右值（multiply()返回值）当左值（matrix&）的错误

本意是想节约空间时间，以免传参时额外复制一份matrix M，那么怎么改才能实现要求呢？

我们不妨利用刚才的右值引用改写上述式子

class matrix{
    matrix(int rows, int lines):
        m_rows(rows), m_lines(lines), m_matrix(new int* [m_lines]){
        for(int i=0; i<m_lines; i++)
            m_matrix[i]=new int [m_rows];
        for(int i=0; i<m_rows; i++)
            for(int j=0; j<m_lines; j++)
                m_matrix[i][j]=0;
    }
    //重载"="
    matrix& operator=(matrix&& M){
        for(int i=0; i<m_lines; i++){
            for(int j=0; j<m_rows; j++){
                m_matrix[i][j]=M.m_matrix[i][j];
            }
        }
        return *this;
    }
    ……
}

int main(){
    matrix M(3,3),A(3.3);
    A=multiply(M,M);
    return 0;
}

运用右值引用成功避免了额外的复制开销，在赋值完成后直接析构掉M。

 

3.如何写出高效的矩阵乘法代码

已知的矩阵乘法算法中，最出名的就是Strassen算法了。它可以将时间复杂度降低到，相比于简单分治和暴力算法的，可以节省大量时间。

原理可以参考CLRS第四章中的说明，简单概括如下：

由于矩阵加减法和矩阵乘法相比计算时间可以忽略不计，即，只要将乘法的计算次数减少，就会降低时间复杂度！

我们先假设矩阵是n*n的方阵，并且n是2的非负整数次幂。

1.普通分治算法（n^3）

在将矩阵分解为四个n/2*n/2的小矩阵后，

，得到A11到A22的计算式：

我们发现有八个乘法，，根据主定理，a=8，b=2，，，所以

可见简单分支并不会改变时间复杂度

2.Strassen分治算法（n^2.807）

原理就是设计了一些矩阵，使得分治后只需计算7个乘法，这样根据主定理，时间复杂度就降低到了。

由于N该方法使用的加法运算较多，如果N<300将不会有明显效果

先创建10个n/2*n/2矩阵S1……S10，令

然后递归计算P1……P7，令

最后可得到结果矩阵C的分块，

我们有18次加法和17个额外存储空间，但只用了7次乘法

附代码实现：

void Strassen_algorithm(int n, const matrix& A, const matrix& B, matrix& R) {
    if (n < 128) {
        R = convention_algorithm(n, n, A, B);
        return;
    }

    matrix S1(n / 2, n / 2), S2(n / 2, n / 2), S3(n / 2, n / 2), S4(n / 2, n / 2), S5(n / 2, n / 2), S6(n / 2, n / 2), S7(n / 2, n / 2), S8(n / 2, n / 2), S9(n / 2, n / 2), S10(n / 2, n / 2), P1(n / 2, n / 2), P2(n / 2, n / 2), P3(n / 2, n / 2), P4(n / 2, n / 2), P5(n / 2, n / 2), P6(n / 2, n / 2), P7(n / 2, n / 2), B11(n / 2, n / 2), B12(n / 2, n / 2), B21(n / 2, n / 2), B22(n / 2, n / 2), A11(n / 2, n / 2), A12(n / 2, n / 2), A21(n / 2, n / 2), A22(n / 2, n / 2), R11(n / 2, n / 2), R12(n / 2, n / 2), R21(n / 2, n / 2), R22(n / 2, n / 2), t1(n / 2, n / 2), t2(n / 2, n / 2);

    int subn = n / 2;
    for (int i = 0; i < subn; i++) {
        for (int j = 0; j < subn; j++) {
            A11.m_matrix[i][j] = A.m_matrix[i][j];
            A12.m_matrix[i][j] = A.m_matrix[i][j + subn];
            A21.m_matrix[i][j] = A.m_matrix[i + subn][j];
            A22.m_matrix[i][j] = A.m_matrix[i + subn][j + subn];
            B11.m_matrix[i][j] = B.m_matrix[i][j];
            B12.m_matrix[i][j] = B.m_matrix[i][j + subn];
            B21.m_matrix[i][j] = B.m_matrix[i + subn][j];
            B22.m_matrix[i][j] = B.m_matrix[i + subn][j + subn];
        }
    }

    minus(B12, B22, S1);
    plus(A11, A12, S2);
    plus(A21, A22, S3);
    minus(B21, B11, S4);
    plus(A11, A22, S5);
    plus(B11, B22, S6);
    minus(A12, A22, S7);
    plus(B21, B22, S8);
    minus(A11, A21, S9);
    plus(B11, B12, S10);

    Strassen_algorithm(subn, A11, S1, P1);
    Strassen_algorithm(subn, S2, B22, P2);
    Strassen_algorithm(subn, S3, B11, P3);
    Strassen_algorithm(subn, A22, S4, P4);
    Strassen_algorithm(subn, S5, S6, P5);
    Strassen_algorithm(subn, S7, S8, P6);
    Strassen_algorithm(subn, S9, S10, P7);

    //R11
    plus(P5, P4, t1);
    minus(t1, P2, t2);
    plus(t2, P6, R11);

    //R12
    plus(P1, P2, R12);

    //R21
    plus(P3, P4, R21);
    
    //R22
    plus(P5, P1, t1);
    minus(t1, P3, t2);
    minus(t2, P7, R22);

    for (int i = 0; i < subn; i++) {
        for (int j = 0; j < subn; j++) {
            R.m_matrix[i][j] = R11.m_matrix[i][j];
            R.m_matrix[i][j + subn] = R12.m_matrix[i][j];
            R.m_matrix[i + subn][j] = R21.m_matrix[i][j];
            R.m_matrix[i + subn][j + subn] = R22.m_matrix[i][j];
        }
    }

    return;
}

很有意思的一点也是很重要的一点，当n小于128时，使用常规算法反而更快。

如果一直到1才开始回退，那么代价实在是太大了，尤其是不断开辟空间释放空间的时间消耗。

在i7上尝试当n=4096时，如果递归到1，十分钟也没跑完，如果递归到64就采用常规算法，那么耗时在30~40秒左右！

然后经过测试，Strassen算法在n为几千时耗时约为暴力算法的1/5左右，且n越大这一比例越小。

同时也测试了不按照内存存储顺序访问的情况，结果连续地址访问比跳跃访问快了约1/10。

3.目前最快的矩阵乘法算法：Coppersmith-Winograd方法（n^2.375）

日后再说吧，在这里放不下，重开一篇详细谈。