AMG工作阶段性总结1

这一段时间，研究了一下的AMG ON GPUS，从一开始的比较乱，应为amgcl用的开源的库还是挺多的，最近把这些知识稍微理了一下，以备下一阶段的学习吧。

首先所有的开始是从http://blog.csdn.net/caiye917015406/article/details/8955928，开始的。

一、用到的相关开源库

   1. Boost库:这个库是amgcl必须的（唯一必须的）。对于boost库，使用时需要自己编译，这里是相关教程：http://blog.csdn.net/caiye917015406/article/details/8957131

   2.OpenCL,如果你需要用GPU运算时，这个是必须的（因为我电脑装的是A卡，所以只能用opencl,当然N卡的话cuda）。在amd官网上下就行了，网上有很多教程。

   3.vexcl:一个方便opencl开发人员使用的库，封装了cl的一些初始化和向量操作。详细的教程：http://blog.csdn.net/caiye917015406/article/details/8976689（这个翻译了一部分）

   4.eigen:这是一个c++的矩阵开源库，包换了许多矩阵操作和数值算法。见：http://blog.csdn.net/caiye917015406/article/category/1417227

   5.viennacl:这是一个支持opencl、openmp和uda的数值开源库，包含了许多实现的数值算法和矩阵操作。这个还没有怎么看，这里有http://viennacl.sourceforge.net/

   当然在上面只有boost库是必须的，其它根据需要来添加。我把所有的都配置好了，以备以后学习用。

二、生产数据

  为了进行测试，我们还需要相关数据，在amgcl里提供了生产测试数据的程序。贴出来吧，我把2D和3D的整理到一块了。

#include <vector>
#include <iostream>
#include <fstream>
#include <cstdlib>

/*
 * Generates problem file for poisson equation in a unit square.
 */
using namespace std;

int gen2D(int N,string filename) {
   // int    n   = argc > 1 ? atoi(argv[1]) : 1024;
	int n = N;
    int    n2  = n * n;
    double h   = 1.0 / (n - 1);
    double h2i = (n - 1) * (n - 1);

    vector<int>    row;
    vector<int>    col;
    vector<double> val;
    vector<double> rhs;

    row.reserve(n2 + 1);
    col.reserve(5 * n2);
    val.reserve(5 * n2);
    rhs.reserve(n2);

    row.push_back(0);

    for (int i = 0, idx = 0; i < n; ++i) {
        double x = i * h;
        for(int j = 0; j < n; ++j, ++idx) {
            double y = j * h;
            if (
                    i == 0 || i == n - 1 ||
                    j == 0 || j == n - 1
               )
            {
                col.push_back(idx);
                val.push_back(1);

                rhs.push_back(0);
            } else {
                col.push_back(idx - n);
                val.push_back(-h2i);

                col.push_back(idx - 1);
                val.push_back(-h2i);

                col.push_back(idx);
                val.push_back(4 * h2i);

                col.push_back(idx + 1);
                val.push_back(-h2i);

                col.push_back(idx + n);
                val.push_back(-h2i);

                rhs.push_back( 2 * (x - x * x + y - y * y) );
            }

            row.push_back(col.size());
        }
    }

    ofstream f(filename, ios::binary);

    f.write((char*)&n2, sizeof(n));
    f.write((char*)row.data(), row.size() * sizeof(row[0]));
    f.write((char*)col.data(), col.size() * sizeof(col[0]));
    f.write((char*)val.data(), val.size() * sizeof(val[0]));
    f.write((char*)rhs.data(), rhs.size() * sizeof(rhs[0]));

    cout << "Wrote \"problem.dat\"" << endl;
	return 0;
}


/*
 * Generates problem file for poisson equation in a unit square.
 */

int gen3D(int N,string filename) {
    int    n   = N;//argc > 1 ? atoi(argv[1]) : 64;
    int    n3  = n * n * n;
    double h2i = (n - 1) * (n - 1);

    vector<int>    row;
    vector<int>    col;
    vector<double> val;
    vector<double> rhs;

    row.reserve(n3 + 1);
    col.reserve(7 * n3);
    val.reserve(7 * n3);
    rhs.reserve(n3);

    row.push_back(0);

    for(int k = 0, idx = 0; k < n; ++k) {
        for(int j = 0; j < n; ++j) {
            for (int i = 0; i < n; ++i, ++idx) {
                if (
                        i == 0 || i == n - 1 ||
                        j == 0 || j == n - 1 ||
                        k == 0 || k == n - 1
                   )
                {
                    col.push_back(idx);
                    val.push_back(1);

                    rhs.push_back(0);
                } else {
                    col.push_back(idx - n * n);
                    val.push_back(-h2i);

                    col.push_back(idx - n);
                    val.push_back(-h2i);

                    col.push_back(idx - 1);
                    val.push_back(-h2i);

                    col.push_back(idx);
                    val.push_back(6 * h2i);

                    col.push_back(idx + 1);
                    val.push_back(-h2i);

                    col.push_back(idx + n);
                    val.push_back(-h2i);

                    col.push_back(idx + n * n);
                    val.push_back(-h2i);

                    rhs.push_back(1);
                }

                row.push_back(col.size());
            }
        }
    }

    ofstream f(filename, ios::binary);

    f.write((char*)&n3, sizeof(n));
    f.write((char*)row.data(), row.size() * sizeof(row[0]));
    f.write((char*)col.data(), col.size() * sizeof(col[0]));
    f.write((char*)val.data(), val.size() * sizeof(val[0]));
    f.write((char*)rhs.data(), rhs.size() * sizeof(rhs[0]));

    cout << "Wrote \"problem.dat\"" << endl;
	return 0;
}

int main()
{
	int n=1024;
	string filename ="E:\\D\\项目\\GPU\\数据\\problem.dat";
	gen2D(n,filename) ;
	return 0;
}

三、读取文件

     在这里生成的是CPS格式的系数矩阵（http://blog.csdn.net/caiye917015406/article/details/8983590）。所以读取的也是相应的矩阵。

#ifndef READ_H
#define READ_H

#include <iostream>
#include <fstream>
#include <vector>
#include <stdexcept>

template <class RHS>
inline int read_problem(const std::string &fname,
        std::vector<int>    &row,
        std::vector<int>    &col,
        std::vector<double> &val,
        RHS &rhs
        )
{
    std::cout << "Reading \"" << fname << "\"..." << std::endl;
    std::ifstream f(fname.c_str(), std::ios::binary);
    if (!f) throw std::inval