1 /* -*- mode: C -*- */
2 /*
3 IGraph library.
4 Copyright (C) 2010-2012 Gabor Csardi <csardi.gabor@gmail.com>
5 334 Harvard st, Cambridge MA, 02139 USA
6
7 This program is free software; you can redistribute it and/or modify
8 it under the terms of the GNU General Public License as published by
9 the Free Software Foundation; either version 2 of the License, or
10 (at your option) any later version.
11
12 This program is distributed in the hope that it will be useful,
13 but WITHOUT ANY WARRANTY; without even the implied warranty of
14 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 GNU General Public License for more details.
16
17 You should have received a copy of the GNU General Public License
18 along with this program; if not, write to the Free Software
19 Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
20 02110-1301 USA
21
22 */
23
24 #include <igraph.h>
25 #include <stdio.h>
26
27 #define DIM 10
28
igraph_print_warning(const char * reason,const char * file,int line,int igraph_errno)29 void igraph_print_warning(const char *reason, const char *file,
30 int line, int igraph_errno) {
31 printf("Warning: %s\n", reason);
32 }
33
main()34 int main() {
35
36 igraph_matrix_t A, B, RHS;
37 int info;
38 int i, j;
39
40 /* Identity matrix, you have to start somewhere */
41
42 igraph_matrix_init(&A, DIM, DIM);
43 igraph_matrix_init(&B, DIM, 1);
44 for (i = 0; i < DIM; i++) {
45 MATRIX(A, i, i) = 1.0;
46 MATRIX(B, i, 0) = i + 1;
47 }
48
49 igraph_matrix_copy(&RHS, &B);
50 igraph_lapack_dgesv(&A, /*ipiv=*/ 0, &RHS, &info);
51
52 if (info != 0) {
53 return 1;
54 }
55 if (!igraph_matrix_all_e(&B, &RHS)) {
56 return 2;
57 }
58
59 igraph_matrix_destroy(&A);
60 igraph_matrix_destroy(&B);
61 igraph_matrix_destroy(&RHS);
62
63 /* Diagonal matrix */
64
65 igraph_matrix_init(&A, DIM, DIM);
66 igraph_matrix_init(&RHS, DIM, 1);
67 for (i = 0; i < DIM; i++) {
68 MATRIX(A, i, i) = i + 1;
69 MATRIX(RHS, i, 0) = i + 1;
70 }
71
72 igraph_lapack_dgesv(&A, /*ipiv=*/ 0, &RHS, &info);
73
74 if (info != 0) {
75 return 3;
76 }
77 for (i = 0; i < DIM; i++) {
78 if (MATRIX(RHS, i, 0) != 1.0) {
79 return 4;
80 }
81 }
82
83 igraph_matrix_destroy(&A);
84 igraph_matrix_destroy(&RHS);
85
86 /* A general matrix */
87
88 igraph_rng_seed(igraph_rng_default(), 42);
89
90 igraph_matrix_init(&A, DIM, DIM);
91 igraph_matrix_init(&B, DIM, 1);
92 igraph_matrix_init(&RHS, DIM, 1);
93 for (i = 0; i < DIM; i++) {
94 int j;
95 MATRIX(B, i, 0) = igraph_rng_get_integer(igraph_rng_default(), 1, 10);
96 for (j = 0; j < DIM; j++) {
97 MATRIX(A, i, j) = igraph_rng_get_integer(igraph_rng_default(), 1, 10);
98 }
99 }
100 igraph_blas_dgemv_array(/*transpose=*/ 0, /*alpha=*/ 1.0, /*a=*/ &A,
101 /*x-*/ &MATRIX(B, 0, 0), /*beta=*/ 0,
102 /*y=*/ &MATRIX(RHS, 0, 0));
103
104 igraph_lapack_dgesv(&A, /*ipiv=*/ 0, &RHS, &info);
105 if (info != 0) {
106 return 5;
107 }
108 for (i = 0; i < DIM; i++) {
109 if (fabs(MATRIX(B, i, 0) - MATRIX(RHS, i, 0)) > 1e-13) {
110 return 6;
111 }
112 }
113
114 igraph_matrix_destroy(&A);
115 igraph_matrix_destroy(&B);
116 igraph_matrix_destroy(&RHS);
117
118 /* A singular matrix */
119
120 igraph_matrix_init(&A, DIM, DIM);
121 igraph_matrix_init(&B, DIM, 1);
122 igraph_matrix_init(&RHS, DIM, 1);
123 for (i = 0; i < DIM; i++) {
124 MATRIX(B, i, 0) = igraph_rng_get_integer(igraph_rng_default(), 1, 10);
125 for (j = 0; j < DIM; j++) {
126 MATRIX(A, i, j) = i == j ? 1 : 0;
127 }
128 }
129 for (i = 0; i < DIM; i++) {
130 MATRIX(A, DIM - 1, i) = MATRIX(A, 0, i);
131 }
132
133 igraph_blas_dgemv_array(/*transpose=*/ 0, /*alpha=*/ 1.0, /*a=*/ &A,
134 /*x-*/ &MATRIX(B, 0, 0), /*beta=*/ 0,
135 /*y=*/ &MATRIX(RHS, 0, 0));
136
137 igraph_set_warning_handler(igraph_print_warning);
138 igraph_lapack_dgesv(&A, /*ipiv=*/ 0, &RHS, &info);
139 if (info != 10) {
140 printf("LAPACK returned info = %d, should have been 10", info);
141 return 7;
142 }
143
144 igraph_matrix_destroy(&A);
145 igraph_matrix_destroy(&B);
146 igraph_matrix_destroy(&RHS);
147
148 return 0;
149 }
150