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
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
|
//------------------------------------------------------------------------------
// GB_is_diagonal: check if A is a diagonal matrix
//------------------------------------------------------------------------------
// SuiteSparse:GraphBLAS, Timothy A. Davis, (c) 2017-2022, All Rights Reserved.
// SPDX-License-Identifier: Apache-2.0
//------------------------------------------------------------------------------
// Returns true if A is a square diagonal matrix, with all diagonal entries
// present. All pending tuples are ignored. Zombies are treated as entries.
#include "GB_mxm.h"
#include "GB_atomics.h"
bool GB_is_diagonal // true if A is diagonal
(
const GrB_Matrix A, // input matrix to examine
GB_Context Context
)
{
//--------------------------------------------------------------------------
// check inputs
//--------------------------------------------------------------------------
ASSERT (A != NULL) ;
ASSERT_MATRIX_OK (A, "A check diag", GB0) ;
ASSERT (!GB_ZOMBIES (A)) ;
ASSERT (GB_JUMBLED_OK (A)) ;
ASSERT (!GB_PENDING (A)) ;
//--------------------------------------------------------------------------
// trivial cases
//--------------------------------------------------------------------------
int64_t n = GB_NROWS (A) ;
int64_t ncols = GB_NCOLS (A) ;
if (n != ncols)
{
// A is rectangular
return (false) ;
}
if (GB_IS_BITMAP (A))
{
// never treat bitmaps as diagonal
return (false) ;
}
if (GB_IS_FULL (A))
{
// A is full, and is diagonal only if 1-by-1, but always return
// false so that GB_AxB_rowscale and GB_AxB_colscale are not used
// by GB_reduce_to_vector.
return (false) ;
}
int64_t anz = GB_nnz (A) ;
int64_t nvec = A->nvec ;
if (n != anz || n != nvec)
{
// A must have exactly n entries in n vectors. A can be sparse or
// hypersparse. If hypersparse, all vectors must be present, so
// Ap has size n+1 whether sparse or hypersparse.
return (false) ;
}
//--------------------------------------------------------------------------
// determine the number of threads to use
//--------------------------------------------------------------------------
// Break the work into lots of tasks so the early-exit can be exploited.
GB_GET_NTHREADS_MAX (nthreads_max, chunk, Context) ;
int nthreads = GB_nthreads (n, chunk, nthreads_max) ;
int ntasks = (nthreads == 1) ? 1 : (256 * nthreads) ;
ntasks = GB_IMIN (ntasks, n) ;
ntasks = GB_IMAX (ntasks, 1) ;
//--------------------------------------------------------------------------
// examine each vector of A
//--------------------------------------------------------------------------
const int64_t *restrict Ap = A->p ;
const int64_t *restrict Ai = A->i ;
int diagonal = true ;
int tid ;
#pragma omp parallel for num_threads(nthreads) schedule(dynamic,1)
for (tid = 0 ; tid < ntasks ; tid++)
{
//----------------------------------------------------------------------
// check for early exit
//----------------------------------------------------------------------
int diag = true ;
{
GB_ATOMIC_READ
diag = diagonal ;
}
if (!diag) continue ;
//----------------------------------------------------------------------
// check if vectors jstart:jend-1 are diagonal
//----------------------------------------------------------------------
int64_t jstart, jend ;
GB_PARTITION (jstart, jend, n, tid, ntasks) ;
for (int64_t j = jstart ; diag && j < jend ; j++)
{
int64_t p = Ap [j] ;
int64_t ajnz = Ap [j+1] - p ;
if (ajnz != 1)
{
// A(:,j) must have exactly one entry
diag = false ;
}
int64_t i = Ai [p] ;
if (i != j)
{
// the single entry must be A(i,i)
diag = false ;
}
}
//----------------------------------------------------------------------
// early exit: tell all other tasks to halt
//----------------------------------------------------------------------
if (!diag)
{
GB_ATOMIC_WRITE
diagonal = false ;
}
}
//--------------------------------------------------------------------------
// return result
//--------------------------------------------------------------------------
return ((bool) diagonal) ;
}
|
