← Back to branch summary

~ubuntu-branches/ubuntu/precise/boinc/precise

« back to all changes in this revision

Viewing changes to samples/nvcuda/cuda_kernel.cu

Tags: 6.12.8+dfsg-1
* New upstream release.
* Simplified debian/rules

Show diffs side-by-side

added added

removed removed

Lines of Context:
samples/nvcuda/cuda_kernel.cu
 
1
// This file is part of BOINC.
 
2
// http://boinc.berkeley.edu
 
3
// Copyright (C) 2008 University of California
 
4
//
 
5
// BOINC is free software; you can redistribute it and/or modify it
 
6
// under the terms of the GNU Lesser General Public License
 
7
// as published by the Free Software Foundation,
 
8
// either version 3 of the License, or (at your option) any later version.
 
9
//
 
10
// BOINC is distributed in the hope that it will be useful,
 
11
// but WITHOUT ANY WARRANTY; without even the implied warranty of
 
12
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
 
13
// See the GNU Lesser General Public License for more details.
 
14
//
 
15
// You should have received a copy of the GNU Lesser General Public License
 
16
// along with BOINC.  If not, see <http://www.gnu.org/licenses/>.
 
17
//
 
18
// This file contains kernel definition for matrix inversion. The external function
 
19
// "invert" serves as an interface between cuda_kernel.cu and cuda.cpp
 
20
//
 
21
// See http://boinc.berkeley.edu/trac/wiki/GPUApp for any compiling issues
 
22
// Contributor: Tuan Le (tuanle86@berkeley.edu)
 
23
 
 
24
// When VERIFY is defined, the sum of squared errors is calculated between the
 
25
// identity matrix and the product A * incerse(A). For debugging...
 
26
//#define VERIFY 1
 
27
#include <stdio.h>
 
28
#include <math.h>
 
29
#include <time.h>
 
30
#include "cuda_config.h"
 
31
 
 
32
__global__ void GEStep1A(REAL * AI, int i, int n2, int lda2) {
 
33
    int k = blockIdx.x * blockDim.x + threadIdx.x;
 
34
    if (k>i && k < n2 && AI[i*lda2+k]!=0) {
 
35
        REAL multiplyer = -AI[i*lda2+k]/AI[i*lda2+i];
 
36
        int n = n2 / 2;
 
37
        for (int j = i+1; j < n; j++) {
 
38
            AI[j*lda2+k] += multiplyer*AI[j*lda2+i];
 
39
        }
 
40
    }
 
41
}
 
42
 
 
43
__global__ void GEStep2(REAL * AI,REAL diag,int i, int n2, int lda2) {
 
44
    int k = blockIdx.x * blockDim.x + threadIdx.x;
 
45
    if (k < n2) {
 
46
        AI[i*lda2+k] /= diag;
 
47
    }
 
48
}
 
49
 
 
50
__global__ void GEStep3(REAL * AI,int i, int n2, int lda2) {
 
51
    int k = blockIdx.x * blockDim.x + threadIdx.x;
 
52
    if (k > i && k < n2) {
 
53
        REAL multiplyer = -AI[i*lda2+k];
 
54
        for (int j = 0; j < i; j++) {
 
55
            AI[j*lda2+k] += multiplyer*AI[j*lda2+i];
 
56
        }
 
57
    }
 
58
}
 
59
 
 
60
/* Helper function for invert. Kernel calls are made in this function */
 
61
void invertge(REAL * AI_d, int lda, int n) {
 
62
    int lda2 = lda * 2;
 
63
    // perform elementary row operations till A in AI becomes identity matrix
 
64
    for (int i = 0; i < n; i++) {
 
65
        GEStep1A<<<(int)ceil((float)(1+(2*n-1)/32)),32>>>(AI_d,i,n*2, lda2);
 
66
        CUDACHECK;
 
67
        cudaThreadSynchronize();
 
68
    }
 
69
 
 
70
    for (int i = n-1; i >= 0; i--) {
 
71
        REAL diag = 1.0;
 
72
        SAFECALL(cudaMemcpy(&diag, &AI_d[i*lda2+i], sizeof(REAL), cudaMemcpyDeviceToHost));
 
73
        GEStep2<<<(int)ceil((float)(1+(n*2-1)/32)),32>>>(AI_d,diag,i,n*2, lda2);
 
74
        CUDACHECK;
 
75
 
 
76
        GEStep3<<<(int)ceil((float)(1+(n*2-1)/32)),32>>>(AI_d,i,n*2, lda2);
 
77
        CUDACHECK;
 
78
        cudaThreadSynchronize();
 
79
        CUDACHECK;
 
80
    }
 
81
}
 
82
 
 
83
/* inverts nxn matrix A and stores result back in A */
 
84
extern void invert(REAL * A, int n) {
 
85
    fprintf(stderr,"starting inversion n = %d ", n);
 
86
    volatile clock_t gputime;
 
87
    gputime=clock();
 
88
 
 
89
    int lda = ((n+15)&~15|16);
 
90
    REAL * AI = (REAL *)malloc(sizeof(REAL)*(n*lda*2));
 
91
    memset(AI,0,sizeof(REAL)*n*lda*2);
 
92
    for (int i = 0; i < n; i++) {
 
93
        memcpy(&AI[lda*i*2], &A[n*i], sizeof(REAL)*n);
 
94
        AI[lda*i*2+n+i] = 1;
 
95
    }
 
96
 
 
97
    REAL * AI_d;
 
98
    SAFECALL(cudaMalloc((void **) &AI_d, sizeof(REAL)*n*lda*2));
 
99
    SAFECALL(cudaMemcpy(AI_d, AI, sizeof(REAL)*n*lda*2, cudaMemcpyHostToDevice));
 
100
 
 
101
    invertge(AI_d, lda, n);
 
102
    SAFECALL(cudaMemcpy(AI, AI_d, sizeof(REAL)*n*lda*2, cudaMemcpyDeviceToHost));
 
103
    cudaFree(AI_d);
 
104
    gputime=clock()-gputime;fprintf(stderr, " %7.1f ms ",gputime/1.e3f);
 
105
    fprintf(stderr, " %7.2f Gflops", 1e-3*(3.0)*n*n*n/3.0/gputime);
 
106
 
 
107
#ifdef VERIFY   
 
108
        // let's verify that
 
109
    REAL error=0.0;
 
110
    // multiply inverse*xcopy, should be Identity matrix
 
111
    for (int k = 0; k < n; k++) {
 
112
        for (int j = 0; j < n; j++) {
 
113
            REAL sum = 0;
 
114
            for (int i = 0; i < n; i++) {
 
115
                sum += AI[j*lda*2+n+i]*A[i*n+k];
 
116
            }
 
117
            if (j!=k) {
 
118
                error += sum * sum;
 
119
            } else {
 
120
                error += (1.0-sum) * (1.0-sum);
 
121
            }
 
122
        }
 
123
    }
 
124
    fprintf(stderr, " %6.2f SSE", error);
 
125
#endif  
 
126
 
 
127
    for (int i = 0; i < n; i++) {
 
128
        memcpy(&A[n*i], &AI[lda*i*2+n], sizeof(REAL)*n);
 
129
    }
 
130
    free(AI);
 
131
    fprintf(stderr," done!\n");
 
132
}
 
 
b'\\ No newline at end of file'