1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
|
/*
* Compute y=A*x
*
* We use
*
*/
#include <stdlib.h>
#include <stdio.h>
#include <csdp/declarations.h>
void matvec(
struct blockmatrix A,
double *x,
double *y)
{
int blk,i,n;
int p;
double *ap;
double scale1;
double scale2;
int inc;
/*
* Work through the blocks one at a time.
*/
p=1;
for (blk=1; blk<=A.nblocks; blk++)
{
switch (A.blocks[blk].blockcategory)
{
case DIAG:
for (i=1; i<=A.blocks[blk].blocksize; i++)
{
y[p]=A.blocks[blk].data.vec[i]*x[p];
p++;
};
break;
case MATRIX:
/*
* Call dgemm to do the matrix multiplication.
*/
n=A.blocks[blk].blocksize;
ap=A.blocks[blk].data.mat;
inc=1;
scale1=1.0;
scale2=0.0;
#ifdef HIDDENSTRLEN
dgemv_("N",&n,&n,&scale1,ap,&n,x+p,&inc,&scale2,y+p,&inc,1);
#else
#ifdef NOUNDERBLAS
#ifdef CAPSBLAS
DGEMV("N",&n,&n,&scale1,ap,&n,x+p,&inc,&scale2,y+p,&inc);
#else
dgemv("N",&n,&n,&scale1,ap,&n,x+p,&inc,&scale2,y+p,&inc);
#endif
#else
#ifdef CAPSBLAS
DGEMV_("N",&n,&n,&scale1,ap,&n,x+p,&inc,&scale2,y+p,&inc);
#else
dgemv_("N",&n,&n,&scale1,ap,&n,x+p,&inc,&scale2,y+p,&inc);
#endif
#endif
#endif
p=p+n;
break;
case PACKEDMATRIX:
default:
printf("matvec illegal block type \n");
exit(206);
};
};
}
|
