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
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
|
//
// lamov.h
//
// Written by Lee Killough 04/19/2012
//
#include "pblas.h"
#include <ctype.h>
extern void xerbla_(const char *, const F_INTG_FCT *, size_t);
void LACPY(const char *UPLO,
const F_INTG_FCT *M,
const F_INTG_FCT *N,
const TYPE *A,
const F_INTG_FCT *LDA,
TYPE *B,
const F_INTG_FCT *LDB);
void LAMOV(const char *UPLO,
const F_INTG_FCT *M,
const F_INTG_FCT *N,
const TYPE *A,
const F_INTG_FCT *LDA,
TYPE *B,
const F_INTG_FCT *LDB)
{
const F_INTG_FCT m = *M;
const F_INTG_FCT n = *N;
const F_INTG_FCT lda = *LDA;
const F_INTG_FCT ldb = *LDB;
if (B + m-1 + ldb*(n-1) < A || A + m-1 + lda*(n-1) < B)
{
LACPY(UPLO, M, N, A, LDA, B, LDB);
}
else if (lda != ldb)
{
TYPE *tmp = malloc(sizeof(*A) * m * n);
if (!tmp)
{
F_INTG_FCT info = -1;
const char func[] = FUNC;
xerbla_(func, &info, sizeof func);
}
else
{
LACPY(UPLO, M, N, A, LDA, tmp, &m);
LACPY(UPLO, M, N, tmp, &m, B, LDB);
free(tmp);
}
}
else
{
F_INTG_FCT i, j;
switch (toupper(*UPLO))
{
case 'U':
if (A > B)
{
for (j=0; j<n; j++)
for (i=0; i<j && i<m; i++)
B[i+ldb*j] = A[i+lda*j];
}
else
{
for (j=n; --j>=0;)
for (i=j<m ? j : m; --i>=0;)
B[i+ldb*j] = A[i+lda*j];
}
break;
case 'L':
if (A > B)
{
for (j=0; j<n; j++)
for (i=j; i<m; i++)
B[i+ldb*j] = A[i+lda*j];
}
else
{
for (j=m<n ? m : n; --j>=0;)
for (i=m; --i>=j;)
B[i+ldb*j] = A[i+lda*j];
}
break;
default:
if (A > B)
{
for (j=0; j<n; j++)
for (i=0; i<m; i++)
B[i+ldb*j] = A[i+lda*j];
}
else
{
for (j=n; --j>=0;)
for (i=m; --i>=0;)
B[i+ldb*j] = A[i+lda*j];
}
break;
}
}
}
|
