| 1619 |
slepc |
1 |
|
|
|
2 |
#include "slepc.h" /*I "slepc.h" I*/
|
|
|
3 |
#include "petscblaslapack.h"
|
|
|
4 |
|
|
|
5 |
/*
|
|
|
6 |
Compute C <- a*A*B + b*C, where
|
|
|
7 |
ldC, the leading dimension of C,
|
|
|
8 |
ldA, the leading dimension of A,
|
|
|
9 |
rA, cA, rows and columns of A,
|
|
|
10 |
At, if true use the transpose of A instead,
|
|
|
11 |
ldB, the leading dimension of B,
|
|
|
12 |
rB, cB, rows and columns of B,
|
|
|
13 |
Bt, if true use the transpose of B instead
|
|
|
14 |
*/
|
|
|
15 |
PetscErrorCode SlepcDenseMatProd(PetscScalar *C, PetscInt _ldC, PetscScalar b,
|
|
|
16 |
PetscScalar a,
|
|
|
17 |
const PetscScalar *A, PetscInt _ldA, PetscInt rA, PetscInt cA, PetscTruth At,
|
|
|
18 |
const PetscScalar *B, PetscInt _ldB, PetscInt rB, PetscInt cB, PetscTruth Bt)
|
|
|
19 |
{
|
|
|
20 |
PetscInt tmp;
|
|
|
21 |
PetscBLASInt m, n, k, ldA = _ldA, ldB = _ldB, ldC = _ldC;
|
|
|
22 |
const char *N = "N", *T = "C", *qA = N, *qB = N;
|
|
|
23 |
|
|
|
24 |
PetscFunctionBegin;
|
|
|
25 |
|
|
|
26 |
/* Transpose if needed */
|
|
|
27 |
if (At == PETSC_TRUE) tmp = rA, rA = cA, cA = tmp, qA = T;
|
|
|
28 |
if (Bt == PETSC_TRUE) tmp = rB, rB = cB, cB = tmp, qB = T;
|
|
|
29 |
|
|
|
30 |
/* Check size */
|
|
|
31 |
if (cA != rB) {
|
|
|
32 |
SETERRQ(1, "Matrix dimensions doesn't match!");
|
|
|
33 |
}
|
|
|
34 |
|
|
|
35 |
/* Do stub */
|
|
|
36 |
m = rA; n = cB; k = cA;
|
|
|
37 |
BLASgemm_(qA, qB, &m, &n, &k, &a, (PetscScalar*)A, &ldA, (PetscScalar*)B,
|
|
|
38 |
&ldB, &b, C, &ldC);
|
|
|
39 |
|
|
|
40 |
PetscFunctionReturn(0);
|
|
|
41 |
}
|
|
|
42 |
|
| 1626 |
slepc |
43 |
/*
|
|
|
44 |
Compute Y[0..cM-1] <- X[0..cX-1] * M, where X and Y are contiguous global
|
|
|
45 |
vectors.
|
|
|
46 |
*/
|
|
|
47 |
PetscErrorCode SlepcUpdateVectorsZ(Vec *Y, Vec *X, PetscInt cX,
|
|
|
48 |
const PetscScalar *M, PetscInt ldM, PetscInt rM, PetscInt cM)
|
|
|
49 |
{
|
|
|
50 |
PetscErrorCode ierr;
|
|
|
51 |
PetscScalar *px, *py;
|
|
|
52 |
PetscInt ldX, ldY;
|
|
|
53 |
|
|
|
54 |
PetscFunctionBegin;
|
|
|
55 |
|
|
|
56 |
if ((cX == 0) || (rM == 0) || (cM == 0)) {
|
|
|
57 |
PetscFunctionReturn(0);
|
|
|
58 |
} else if ((Y + cM <= X) || (X + cX <= Y)) {
|
|
|
59 |
/* If Y[0..cM-1] and X[0..cX-1] are not overlapped... */
|
|
|
60 |
|
|
|
61 |
/* Get the dense matrices and dimensions associated to Y and X */
|
|
|
62 |
ierr = VecGetLocalSize(X[0], &ldX); CHKERRQ(ierr);
|
|
|
63 |
ierr = VecGetLocalSize(Y[0], &ldY); CHKERRQ(ierr);
|
|
|
64 |
if ((ldX != ldY) || (cX != rM)) {
|
|
|
65 |
SETERRQ(1, "Matrix dimensions doesn't match!");
|
|
|
66 |
}
|
|
|
67 |
ierr = VecGetArray(X[0], &px);CHKERRQ(ierr);
|
|
|
68 |
ierr = VecGetArray(Y[0], &py);CHKERRQ(ierr);
|
|
|
69 |
|
|
|
70 |
/* Do operation */
|
|
|
71 |
ierr = SlepcDenseMatProd(py, ldX, 0.0, 1.0, px, ldX, ldX, cX, PETSC_FALSE,
|
|
|
72 |
M, ldM, rM, cM, PETSC_FALSE); CHKERRQ(ierr);
|
|
|
73 |
|
|
|
74 |
|
|
|
75 |
ierr = VecRestoreArray(X[0], &px);CHKERRQ(ierr);
|
|
|
76 |
ierr = VecRestoreArray(Y[0], &py);CHKERRQ(ierr);
|
|
|
77 |
|
|
|
78 |
} else if (Y > X) {
|
|
|
79 |
/* If not, call to SlepcUpdateVectors */
|
|
|
80 |
ierr = SlepcUpdateVectors(cX, X, Y-X, Y-X+cM, M, ldM, PETSC_FALSE);
|
|
|
81 |
CHKERRQ(ierr);
|
|
|
82 |
} else {
|
|
|
83 |
SETERRQ(1, "I don't support this case!");
|
|
|
84 |
}
|
|
|
85 |
|
|
|
86 |
PetscFunctionReturn(0);
|
|
|
87 |
}
|
|
|
88 |
|